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Abstract 


The  sum-difference  surface  current  formulation  is  introduced  to  treat  electromagnetic 
boundary  value  problems  when  anisotropic  impedances  are  specified  on  both  sides  of  a 
surface.  It  can  also  be  applied  to  impedance  coated  bodies.  This  formulation  preserves  the 
duality  nature  of  Maxwell  equations  and  carries  it  over  into  the  algebraic  form  of  the 
integrodifferential  operators  in  the  equations  for  surface  currents.  Since  a  90°  rotation  is 
equivalent  to  undergoing  a  duality  transform  for  an  incident  plane  wave,  this  particular 
symmetry  in  the  algebraic  form  of  the  operators  leads  to  sufficient  conditions  under  which 
the  on-axis  backscattering  of  an  anisotropic  impedance  coated  scatterer  having  a  90° 
rotational  symmetry  is  eliminated.  The  sum-difference  formulation  is  utilized  for  solving  the 
problem  of  electromagnetic  scattering  from  an  anisotropically  impedance  coated  tubular 
cylinder  of  finite  length.  The  solution  has  been  coded  in  FORTRAN  and  tested.  Some 
interesting  results  are  presented  and  discussed. 
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I.  INTRODUCTION 


Sometime  ago  the  question  was  raised:  "For  electromagnetic  boundary  value 
problems  with  specified  surface  impedances,  how  can  one  go  from  a  non-perfectly 
conducting  surface  on  which  both  the  electric  and  the  magnetic  equivalent  surface  currents 
are  to  be  found,  to  a  perfectly  conducting  surface  on  which  the  number  of  unknowns  is 
halved  [1]?"  The  answer  to  this  question  turns  out  to  be  one  of  algebra.  It  is  well  known  that 
the  impedance  specified  on  the  surface  of  a  body  separates  its  interior  completely  from  its 
exterior.  Therefore  an  impedance  coated  body  can  always  be  considered  as  a  hollow  volume 
enclosed  by  an  infinitesimally  thin  shell  with  surface  impedances  specified  both  on  the  inside 
and  the  outside  of  the  shell.  The  inside  and  the  outside  of  the  body  can  be  considered  as 
constituted  of  the  same  medium  and  the  impressed  electromagnetic  excitation  can  be  treated 
as  continuous  across  the  shell.  On  the  outside  surface,  there  are  the  equivalent  total  electric 

current  K +  and  total  magnetic  current  L  * ;  on  the  inside  surface,  there  are  K  ’  and  L  .  For 
an  exterior  problem,  only  K +  and  L +  need  to  be  found;  for  an  interior  problem,  only  K  ‘  and 

L  are  necessary.  A  single  formulation  for  solving  both  types  of  problems  would  appear  to 

require  finding  all  inside  and  outside  currents  therefore  doubling  the  amount  of  work,  but  it 
turns  out  not  to  be  the  case  because  some  of  the  currents  are  linear  combinations  of  others. 
Furthermore,  this  formulation  holds  the  key  to  answering  the  question  posed  above. 

Since  the  shell  is  infinitesimally  thin,  from  Maxwell  equations  the  radiation  to  the 
outside  and  to  the  inside  of  the  shell  can  both  be  given  in  terms  of  integrodifferential 

operators  on  the  sum  currents  K  =  K+  +  K~  and  L  =  L +  +  L  ” .  Note  that  the  outside  currents 

will  not  contribute  to  the  radiation  in  the  interior  while  the  inside  currents  will  not  contribute 
to  the  radiation  to  the  exterior  of  the  shell.  For  simplicity  in  the  description,  we  consider  the 
exterior  problem  of  electromagnetic  scattering.  By  definition,  the  radiation  is  the  difference 
between  the  total  field  and  the  incident  field.  On  the  surfaces  of  the  shell,  this  definition  links 

the  incident  E  and  H  fields  and  the  difference  currents  K-K  and  L  *  -L  to  the  sum 
currents. Tt  is  therefore  natural  to  treat  the  difference  currents  and  the  sum  currents  as  the  four 
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unknowns  to  be  solved  instead  of  the  inside  and  the  outside  currents. 

The  surface  impedance  on  the  outside  surface  of  the  shell  normalized  to  the  medium 
is  denoted  by  Z+  and  that  on  the  inside  surface  is  Z\  They  can  be  tensors  if  the  impedances 
are  anisotropic  and  may  vary  from  point  to  point.  By  forming  the  sum  impedance 
Z  =  (Z +  +  Z ' )/ 2  and  the  difference  impedance  A  =  (Z +  -  Z  ~  )/2 ,  the  impedance  boundary 
conditions  provide  a  set  of  relations  between  the  difference  and  the  sum  currents.  It  turns  out 
that  the  rank  of  Z  determines  how  the  unknown  surface  currents  are  solved.  If  Z  is  invertible, 
then  the  difference  currents  are  linear  combinations  of  the  sum  currents  so  that  only  the 
integrodifferential  equations  of  the  sum  currents  have  to  be  solved.  There  are  only  two 
unknowns  to  be  solved  for  both  the  exterior  and  the  interior  problems  under  this  sum- 
difference  current  formulation.  If  Z  is  rank  0,  then  Z  =  0;  the  impedance  boundary  condition 
requires  that  L  be  proportional  to  a  90°  rotation  of  A.K.  The  difference  electric  current  is 

obtained  from  K  after  the  integrodifferential  equation  on  K  is  solved.  This  situation  includes 

the  case  where  Z +  =  Z  *  =  0  when  the  surface  is  perfectly  conducting,  thus  the  result  answers 
the  question  about  the  transition  of  the  equations  from  a  problem  of  two  variables  to  one 
which  has  only  a  single  variable. 

Instead  of  dealing  with  an  impedance  coated  body,  this  thesis  presents  the  sum- 
difference  currents  formulation  of  electromagnetic  boundary  value  problems  for  the 
scattering  of  an  infinitesimally  thin  surface  for  which  both  the  inside  and  the  outside  currents 
are  true  unknowns  to  be  found.  Extension  of  this  formulation  to  impedance  coated  bodies  is 
then  discussed. 

This  formulation  preserves  the  duality  nature  of  Maxwell  equations  and  carries  it  over 
into  an  explicit  specific  algebraic  form  of  the  integrodifferential  operators  in  the  equations 
for  the  sum  currents.  Since,  for  a  plane  wave,  a  90°  rotation  is  equivalent  to  undergoing  a 
duality  transform,  this  explicit  symmetry  in  the  algebraic  form  of  the  operators  enables  us  to 
deduce  sufficient  conditions  for  eliminating  the  on-axis  backscattering  of  an  anisotropic 
impedance  coated  scatterer  having  a  90°  rotational  symmetry. 

The  sum-difference  currents  formulation  is  utilized  for  solving  the  problem  of 
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electromagnetic  scattering  from  an  anisotropically  impedance  coated  tubular  cylinder  of 
finite  length.  The  solution  has  been  coded  in  FORTRAN  and  tested.  Some  interesting  results 
are  presented  and  discussed. 

In  this  thesis,  the  time  dependence  e is  used.  E  represents  the  electric  field 
intensity  divided  by  the  intrinsic  impedance  of  the  medium,  y/p7e ;  therefore  it  has  the  same 
unit  as  IT  in  amperes  per  meter.  So  are  the  electric  and  magnetic  equivalent  surface  currents  K 

and  L. 


3 


II.  THE  SUM-DIFFERENCE  SURFACE  CURRENT  FORMULATION  OF 
ELECTROMAGNETIC  BOUNDARY-VALUE  PROBLEMS 

A.  STRATTON-CHU  FIELD  FORMULATION  AND  RADIATION 

On  an  orientable,  piecewise  regular  surface,  whether  open  or  closed,  having  a  surface 
electric  current  density  K  and  a  surface  magnetic  current  density  L,  we  define  the  Stratton- 
Chu  E-field  formula  as: 

EsJr,S,K,L)  =  ^  f  K(ro)G(T-ro)da0-  -j-V  f  K(rc)-V0G(?-rc)da0  . 

4tc  J  s  4tc  **  s  ^ 

-  ^  X  ffivar-rjda. 

ik\r-T0\ 

- ,  then  the 

r~r0\ 

Stratton-Chu  H-field  formula  can  be  defined  as: 

HsJr,S,K,L)  =  Es_c(r,S,L,-K) 

-  ^  X  [g(V<xr-rJda.  *  U  fuocur-fjda.  (2.2) 

Note  that  if  S  is  a  closed  surface  and  K  and  L  are  the  actual  total  equivalent  surface  currents 
on  S,  then  Es_c  and  Hs_c  are  respectively  the  E  and  H  fields  at  r  due  to  all  sources  inside 
S  if  7  is  located  outside  S  and  vice  versa.  This  is  a  direct  consequence  of  Maxwell’s 
equations  [2]  and  under  this  circumstance,  the  Stratton-Chu  formulae  are  equivalent  to 
Maxwell’s  equations.  On  the  other  hand,  unlike  the  Poynting  theorem,  the  Stratton-Chu  field 
formulae  over  an  open  surface  S  are  but  integrodifferential  operators  on  the  tangential  vector 


where  r  is  a  point  which  is  not  on  S,  k  =  and  G{r-rg)  = 
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fields  K  and  L  over  S  without  any  special  physical  meaning  attached. 

To  introduce  equivalent  surface  currents  on  S,  the  direction  of  the  unit  normal  vector  n 
on  the  surface  has  to  be  chosen.  Adopting  the  terminology  for  a  closed  surface,  we  can  assign 
one  side  of  any  orientable  surface  S  as  the  "outside  surface"  S+,  albeit  somewhat  arbitrarily 
if  S  is  not  closed.  The  outward  normal  n  *  is  the  unit  normal  pointing  out  of  this  side  of  S 
and  for  simplicity,  we  call  n  +  the  normal  of  S  and  denote  it  by  n.  The  other  side  of  the 
surface  S  is  the  "inside  surface"  S'.  At  any  point  r  on  S,  the  inward  normal  n  '  is  -n.  As  a 

convention,  the  fields  and  surface  currents  on  S  +  and  S "  always  carry  the  corresponding 
superscripts  (Figure  2-1). 


Figure  2-1.  Outside  and  Inside  Surfaces,  Normals  and  the  Equivalent  Currents. 

Each  of  the  total  surface  currents  K*  and  L±  on  S±  consists  of  two  parts:  the  incident 
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current  (with  the  additional  superscript  "me")  and  the  scattering  current  (with  the  additional 
superscript  "sc")  corresponding  to  the  incident  field  and  the  scattered  field  on  the  particular 
side  of  the  surface  S: 


K*"  =  ^  x  H*  =  n*  x  Hmc  +  n*  x  H 


±,SC 


=  K*  +  K' 


(2-3) 


L  —  7— <  —  ^  -4-  7^  IflC  /\  4-  7-!  H,  SC  /s  -f- 

=  E  x  =  E  x  n~  +  E  x  n~ 


L  ±,  me  7  : 

+  L 


(2-4) 


Note  that  5  is  infinitesimally  thin,  hence  H +’mc  =  H  ~'inc 


fr  inc  j  f*+,inc  ^-Jnc  inc 

=  H  and  E  =E  =  E  on 


S  so  that  K  ,inc  =  - K  ’inc  and  L*'inc  =  -  L  'mc .  Therefore  the  sum  currents  K  and  L  on  S 
defined  below  are  also  the  corresponding  sums  of  the  scattering  currents  only: 


K  =  K+  +  K  =  K+’sc  +  K~'sc 


L  =  L  +  L  =  L  +  L 


*  ~,SC 


(2-5) 


Since  the  Stratton-Chu  field  formulae  are  linear  operators  on  the  surface  currents,  the 
radiated  fields  from  surface  currents  on  S  are  determined  by  the  sum  currents  only: 

Es\f)  =  Es_c{r,S\K\C)  +  EsJf,S-,k'X ')  =  EsJr,S,K,L)  _ 

(2-6) 

Hx(r)  =  HsJf,S,K,L)  =  Es_c(r,S,L,-K ) 

B.  CONDITION  ON  THE  CURRENTS  IMPOSED  BY  MAXWELL’S 
EQUATIONS 

As  7  approaches  7±  on  S*,  the  tangential  components  (denoted  by  the  subscript 
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"tan")  of  Eq.  (2-6)  provide  four  equations  relating  the  incident  fields  and  the  total  currents 
on  both  sides  of  S  through  the  fact  that  the  incident  field  is  the  difference  between  the  total 
and  the  scattered  field: 


C  =  Cn  -  Kc,  JT,SX£) 
C  =  C  -  KcJT’SXl) 
CC  =  -  Hs.cm(r\S,K,L) 

Cc  =  C  -  Hs.c^(r,S,K,L) 


(2-7) 

(2-8) 

(2-9) 

(2-10) 


These  four  equations  are  not  independent  of  each  other.  Because 
n±  x  (£(?*)  x  n*)  =  =  n*  x  L* 

n ±  x  (^(r4)  x  n *)  =  ^(r*)  =  -n±  x  f ± 


(2-11) 


the  difference  between  Eqs.  (2-7)  and  (2-8)  trivially  confirms  the  definition  of  the  sum 
equivalent  magnetic  current  while  the  difference  between  Eqs.  (2-9)  and  (2-10)  confirms  the 
definition  of  the  sum  equivalent  electric  current  as  both  can  be  deduced  directly  from 
Maxwell’s  equations.  We  choose  to  use  the  sum  of  Eqs.  (2-7)  and  (2-8)  and  that  of  Eqs.  (2-9) 
and  (2-10)  as  the  two  independent  linear  combinations  out  of  Eqs.  (2-7)  through  (2-10)  to 
link  the  incident  fields  to  the  total  surface  currents  on  S*  as  dictated  by  Maxwell’s  equations: 
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2E j  =  n  x  (L+  -  L)  -  { EsJr+,S,K,L] )  +  Esjr,S,K,L)} 

=  n  x  (L+  -  L")  +  -  NL 

2HZ  =  -Ax(K*  -  K~)  -  {^_c(f+,S,2?,L)  +  tf,_c(r,S,£L)} 
=  -  n  x  (K+  -  K~)  +  NK  +  ML 


(2-12) 


where  M  and  N  are  linear  integrodifferential  operators  on  the  tangential  vector  fields  K  and 
L  over  S. 

Under  any  orthonormal  coordinates  ( u ,  v)  over  S  having  u,  v  as  the  unit  basis  vectors 
and  with  n  =  u  x  v,  a  tangential  vector  field  A  over  S  can  be  written  in  matrix  form  as: 

is  one  of  the  Pauli  spin 


A  = 

.  Then  n  x  A  = 

=  -io~A  where  o,  = 

0  -i 

.  v 

2  2 

i  0 

matrices.  Note  that  o22  =  1.  Using  these  matrix  notations,  we  can  rewrite  Eq.  (2-12)  in  the 
following  form: 


0  ia2 

r 

-K~ 

M  -N 

K 

i ?inc 

^tan 

- 

-  2 

-ia2  0 

r 

-L 

£ 

L 

H.inc 

L  J 

L 

tan 

(2-13) 


C.  IMPEDANCE  BOUNDARY  CONDITION 

Maxwell’s  equations  alone  cannot  determine  the  electromagnetic  fields  completely. 
If  S  is  an  open  surface,  appropriate  boundary  condition  which  the  fields  satisfy  on  S  must  be 
specified.  It  is  usually  given  in  terms  of  the  impedance  boundary  condition,  a  linear  relation 

among  the  tangential  components  of  the  total  E  and  the  total  H  fields  on  S.  If  S  is  a  closed 
surface,  there  are  two  possibilities:  One  is  to  specify  the  electric  and  magnetic  properties  of 
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the  volume  within  S  and  require  the  fields  to  satisfy  regularity  conditions  within  S  and  be 
linked  to  the  fields  outside  through  boundary  conditions  across  S;  another  is  to  specify  the 
impedance  boundary  condition  on  S  +  for  an  exterior  problem  or  on' 5  for  an  interior 
problem.  Note  that  an  impedance  boundary  condition  over  a  closed  surface  S  completely 
separates  the  exterior  from  the  interior  of  S.  Therefore,  the  surface  impedance  on  S'  can  be 
arbitrary  for  an  exterior  problem  while  that  on  S  +  can  be  arbitrary  for  an  interior  problem. 
In  this  thesis,  normalized  surface  impedances  Z±  are  assumed  to  be  specified  on  S  ±  whether 
S  is  an  open  or  a  closed  surface.  Note  that  an  open  surface  S  can  be  considered  as  bounded 
within  the  closed  surface  formed  by  joining  S+  and  S  \ 

The  impedance  boundary  conditions  on  S±  are  defined  by: 

A*  x  (£*  x  n±)  =  Z±  (rf*  x  H*)  (2-14) 


or  equivalently,  in  terms  of  the  total  surface  currents: 

n±  x  L±=  Z*  K 4  (2-15) 

With  the  matrix  notations  for  tangential  vector  fields  over  S  in  the  orthonormal  coordinate 
system  («,  v)  introduced  before,  we  can  consider  Z±  as  2x2  matrices  and  rewrite  Eq.  (2-15) 
in  the  form: 


+  io2  L±  =  Z±K±=  ^Z±[K±  (K+  -  K')] 


(2-16) 


which  can  readily  be  recast  into  a  relation  among  sum  and  difference  currents: 


- 1 

1 

> 

1 

Q 

to 

K  + 

-  K 

z 

0 

K 

z  0 

V 

-  V 

-A 

1 

Q 

.to 

L 

(2-17) 
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with 


z  =  1  (Z+  +  Z-)  (2-18) 

and 

A  =  1  (Z*  -  Z-)  (2-19) 

Eqs.  (2-13)  and  (2-17)  are  a  set  of  four  two-dimensional  vector  equations  to  be  solved  for  the 
sum  and  difference  equivalent  electric  and  magnetic  surface  current  densities  on  S. 

I).  ALGEBRA  OF  THE  SUM-DIFFERENCE  CURRENT  EQUATIONS 

The  existence  and  uniqueness  of  a  solution  to  either  the  exterior  or  the  interior 
problem  specified  in  terms  of  the -impedance  boundary  condition  have  been  well  established 
[3],  Here  we  want  to  investigate  how  such  a  solution  can  be  obtained  from  Eqs.  (2-13)  and 
(2-17).  Clearly  Eq.  (2-17)  defines  uniquely  the  algebraic  relationship  between  the  difference 
and  the  sum  currents  if  Z  is  invertible.  For  example,  the  difference  currents  can  be  expressed 
in  terms  of  the  sum  currents: 

-i 

K*  -  K~  0  i°2]  \k 

R  (2-20) 

V-L  \  [-io2  °j  [L 

where 

Z- AZ-1  A  -i  AZ'1  o2 

R  =  (2-21) 

-io2Z  *A  o2Z  *a2 
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An  equation  for  the  sum  currents  is  obtained  by  substituting  Eq.  (2-20)  into  Eq.  (2-13): 


r 

pinC 

^tan 

M 

-N 

K 

+  R 

K 

=  2 

N 

M 

L 

L 

TJ  inC 

(2-22) 


which  can  be  solved  for  K  and  L.  Eq.  (2-20)  in  turn  enables  us  to  compute  the  difference 

currents  algebraically  then  split  the  sum  and  the  difference  currents  into  total  currents  on  S*. 

If  Z  is  not  invertible,  then  the  situation  is  more  complicated.  Z  can  either  be  of  rank 
0  when  Z  =  0  or  rank  1  when  det  Z  =  0  but  Z  *■  0.  Eqs.  (2-13)  and  (2-17)  can  be  combined 

into  an  equation  for  the  sum  currents  K  and  L: 


1 

0 


i  A  o2 
—iZo  2 


K 

N  M 

L 

Z 

-A 


K 

L 


iAo2 

-iZ 


tan 


H. 


tan 


(2-23) 


When  Z  =  0,  Eq.  (2-17)  gives  the  null  relations  L  +  -  L~  =  io2  A  (K*  -  K~)  and 
L  =  i  o2  A  K  hence  L~  =  i  o2  A  K±.  Eq.  (2-23)  becomes  one  for  K  only: 


M  -N 

1 

Einc 

■^tan 

[l  i  A  a J 

K  =  2  fl  i  A  o ,1 

L  2J 

N  M 

io2  A 

1  2J 

Hinc 

tan 

Eq.  (2-13)  has  to  be  used  to  find  the  difference  electric  current: 

K+  -  K~  =  ioz[N  +  iMo2 A]  K  -  2 ia2H^  (2-25) 


Therefore, 
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(2-26) 


P  =  -i{l  ±  io2[W  *  iMo2 A]  }*  t  io2fi“ 

Since  the  last  term  in  Eq.  (2-26)  is  A?±,wc , 

=  I{1  ±  /a2[jV  +  iMo2A]  }£  (2-27) 


L +  and  L  can  be  obtained  algebraically  by  multiplying  io2  A  to  K +  and  K  respectively. 
On  the  other  hand,  by  Eq.  (2-13), 

L±,sc  =  ^ia2{A  t  [M  -  iWo2A]  }K  (2-28) 

Note  that  the  Z  =  0  case  includes  the  special  situation  Z+  =  Z'  =  Z=  A  =  0  when 
both  sides  of  S  are  perfectly  conducting.  Under  this  circumstance  L  =  LT  =  0  and  the 
operator  N  is  never  involved. 

When  Z  is  rank  1 ,  Z  *  0  but  det  Z  =  0.  The  right  hand  side  of  Eq.  (2- 1 7)  provides  one 
linear  relation  between  the  components  of  L  -  io2AK  which  can  be  used  to  reduce  the  four 

components  of  K  and  L  as  the  unknown  quantities  in  Eq.  (2-23)  to  three  so  that  the 
remaining  three  components  can  be  solved.  The  left  hand  side  of  Eq.  (2-17)  assures  that  the 
same  linear  relationship  between  the  components  of  L-io2AK  exists  _between 

corresponding  components  of  the  difference  currents.  Eq.  (2-13)  again  has  to  be  invoked  to 
compute  three  other  linearly  independent  combinations  of  the  components  of  the  difference 
currents  from  the  sum  currents. 
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E.  CONSIDERATIONS  FOR  A  CLOSED  SURFACE 


When  S  is  a  closed  surface,  the  choice  of  Z'  can  be  arbitrary  for  an  exterior  problem 
such  as  scattering  while  the  choice  of  Z+  is  arbitrary  for  an  interior  problem.  It  is  desirable 
to  choose  Z"  =  -Z+  so  that  Z  =  0  and  A  =  Z+  =  -Z'.  Then  we  have  L  =  io2AK  and  only  K 

has  to  be  computed.  With  such  a  choice,  for  an  exterior  problem, 


K*’sc  =  1(1  -  o2Mo2Z+  +  io2N)K 


(2-29) 


iNo2Z  +  -  M)K 


(2-30) 


and  for  an  interior  problem: 

r-“  »  1(1  -  o2Mo2Z-  -  ia2N)K  (2-31) 


L',SC  =  -  Z-  +  i N  o2Z')K  (2-32) 
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III.  A  THEOREM  OF  ANISOTROPIC  ABSORBERS 


A.  AXIAL  RADIATION  FROM  A  SURFACE  OF  90°  ROTATIONAL 
SYMMETRY 


The  ry-plane  cross  section  of  a  surface  S  having  a  90  °  rotational  symmetry  around 
the  z-axis  is  shown  in  Figure  3-1.  Because  of  this  symmetry,  S  can  be  decomposed  into  four 
non-overlapping  congruent  pieces  S„  S2,  S3,  S4  so  that  a  90  °  rotation  will  bring  5,  into  S,+!. 


Figure  3-1.  Cross  Section  of  a  Surface  of  90-Degree  Rotational  Symmetry. 


(These  subscripts  are  considered  as  equal  under  modulo  4.)  Therefore  each  piece  S.  of  S  can 
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dr.  dr. 

be  parametrized  in  the  same  orthonormal  coordinates  ( u ,  v),  with  u  =  — ,  v  =  —  the 

du  dv 


orthonormal  basis  vectors  on  St,  as  follows: 

r,  =  (xo(u,v),yo(u,v),zo(u,v)) 
r f2  =  (-yo(u,v),xo(u,v),zo(u,v) ) 

(3-D 

r3  =  (-xg(u,v),-yo(u,v),z0(u,v)) 
rA  =  (y0(u,v),  -xo(u,v),zo(u,v)) 

where  r.  e  Sr  As  an  example  of  a  possible  choice,  u  =  constant  and  v  =  constant  can  be  the 
lines  of  curvature  of  Sj. 

In  terms  of  the  coordinates  (u,  v),  the  sum  surface  current  distributions  on  S,  are: 

K(F)  =  K.(u,v) 

(3-2) 

L(f)  =  L.(«,v) 


Since  for  r  »  ro. 


G(r)  - 


e 


I 


kr 


(3-3) 


VG(F)  «  ikrG  «  -VtfG(r) 


(3-4) 


the  radiation  from  such  current  distributions  at  a  distance  r  »  max  |  r,  |  along  the  positive 
z-axis  is,  from  Eq.  (2-6): 
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Esc(r)  =  EfJr,S,K,L ) 
ik 

4tc  rJs 
ik 

4  nr 


fs{*WV  +L/01  +MKy(0 -Lx(f0)]}eik^^dao 

f  /s|^E[^v)+L.(m,v)] 

-  st  [^v)-l>,v)]L^^-^^^2  dudv 


(3-5) 


Note  that  this  approximation  is  independent  of  the  wavelength;  it  is  applicable  in  regions 
closer  to  S  than  the  usual  Fresnel  zone. 

B.  CONDITION  FOR  VANISHING  ON-AXIS  BACKSCATTERING 

Consider  two  situations  when  a  linearly  polarized  plane  wave  of  unit  strength  is 
incident  on  S  along  the  z-axis  from  the  positive  direction;  Situation  1,  identified  with  the 
superscript  (1)  has  the  wave  polarized  in  the  x-direction  while  Situation  2,  identified  with  the 
superscript  (2),  has  the  wave  polarized  in  the  y-direction.  The  incident  fields  are  respectively: 


EmcXl)  =  xe~ikz 

gte’W  =  _  ^  e-ikz 


(3-6) 


EincX 2)  =  y  e~ikz 

HincA2)  =  f  e  - ikz 


Note  that,  as  seen  from  the  positive  z-axis,  the  incident  wave  in  Situation  2  is  that  of 
Situation  1  rotated  by  90°  counterclockwise.  Furthermore,  Situation  2  can  be  obtained  from 

Situation  1  through  the  duality  transformation  Emc  -  H  mc,  Hinc  —  E  mc .  Therefore,  for 

a  plane  wave,  a  90°  rotation  is  equivalent  to  undergoing  the  duality  transform. 

Because  of  the  rotational  symmetry  of  S,  the  currents  excited  on  5,  under  Situation 
1  must  appear  on  SM  under  Situation  2.  Therefore: 
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(3-8) 


K^[x(u,v)  =  -  ^(«,v) 

K™hy(u,v)  =  K^(u,v) 

K«\ju,v)  =  K™(u,v) 

L^\  x(u,v)  =  -  L-  ly\u,v) 

L^\y(u,v)  =  L-'x\u,v)  (3-9) 

Lffu(u,v)  =  L.  lz\u,v) 


Assume  that  Z  on  S  is  invertible,  the  sum  currents  on  S  are  determined  by  Eq.  (2-22).  The 
tangential  components  of  the  incident  fields  which  appear  on  the  right-hand-side  of  that 
equation  under  these  two  situations  are: 


inc,  ( 1 ) 
^tan 


H 


inc,(  1) 


tan 


U'X 

A  A 

v-x 


e-ikz 


- u-y 

A  /S 

-vy 


(3-10) 


77. 


line,  (2) 
'tan 

inc,  (2) 
tan 


A  A 

wy 


v-y 

11 

i 

0 

-7 

ffinc,(  1) 
^  tan 

u-x 

7 

0 

rjMC,(  1) 

_  tan 

V'X 


(3-11) 


where  7  is  the  2x2  identity  matrix.  Therefore, 
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(3-12) 


M  -N 

KW 

N  M 

L(1) 

=  2 


£?<nc,( 1) 
^tan 


Trine,  (1) 
^tan 


M 

N 


-N 

M 


+  R 


=  2 


ffinc,(2) 

^tan 


inc,(2) 

tan 


pineal) 

^tan 

* 

rjinc,(l) 
17  f 
tan 

(3-13) 


Since 


1 

O 

M  -N 

I  0 

N  M 

M  -N 

*-< 

1 

o 

N  M 

1 - 

o 

1 _ 

•  (3-14) 


it  follows  that  if 


0 

I 


0 

/ 


-/ 

0 


(3-15) 


we  can  multiply 


0 

I 


-I 

0 


to  Eq.  (3-12)  to  get: 


M  -N 

o 

1 - 

KW 

+  R 

0 

'] 

Km 

1! 

N  M 

U 

0  ' 

L(1) 

I 

0 

L 

1 

pinc,(\) 

tan 

- 

frinc  A 1) 
^tan 

(3-16) 


Therefore  the  excited  surface  currents  in  these  two  situations  are  related  by: 


o  ’ 

i 

>--< , 

!*U 

'W' 

_ 1 

-L(1) 

r® . 

L 1  °. 

lf». 

£d) 

(3-17) 
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Combining  this  result  with  Eqs.  (3-8)  and  (3-9),  we  have: 


,>,v)  =  -  L^]\x(u, v)  =  -  K™(u,v) 
K™  y(u,v)  =  -  L-l\  y(u,v)  =  K£\u,v) 
L™hx(u,v)  =  K;l\x(u,v)  =  -  L$(u,v) 
Lil i,y(«.v)  =  ^i,/«,v)  =  L-'x\u,v) 

so  that 

E  &£(«,*)  +  L,^  («,v)]  =  0 

/  =  1 

and 


E  [<?(«.  v)  -  L^fav)]  =  0 


(  =  i 


Hence,  along  the  positive  z-axis,  from  Eq  (3-5), 
ik 


E  (r)  = 


4-Tcr 


*E  [*£W)  +  L«\u,v)] 


l~  1 

4  r 

+  *E 

i— 1 


0«.v)  ■  LS(u’v) 


ik\J(z-zf+x20+yl 


dudv 


=  0 


and  the  backscattering  from  S  along  the  positive  z-direction  must  vanish  if  Eq. 
satisfied. 


(3-18) 


(3-19) 


(3-20) 


(3-21) 


(3-15)  is 
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C.  IMPEDANCE  MATRICES  FOR  ZERO  ON-AXIS  BACKS  CATTERING 

Z-AZ~!A  -i  AZ-1o, 


It  can  be  verified  that  the  matrix 


-io2  Z'1  A  o2  Z  ‘o2 


commutes  with 


0  -I 
I  0 


if  and  only  if: 


o2Z‘,A  =  -  A  Z"1o2 


(3-22) 


Z  -  a2  Z^Oj  =  AZ_1A 


(3-23) 


where  both  Z  and  A  are  2  x  2  matrices.  Under  the  assumption  that  the  inverse  of  Z  exists, 
we  analyze  Eqs.  (3-22)  and  (3-23)  as  follows: 

Because  of  the  identity: 


-i 


detZ 


a  2  Z 


a. 


(3-24) 


and,  by  multiplying  o2  to  both  sides  of  Eq.  (3-22): 

Z’1  A  =  -  o2  A  Z~l  o2  (3-25) 


Eq.  (3-23)  can  now  be  transformed  to 


-  A  (  o2  A  Z"1  o2  )  +  o2  Z  1  o2 

-  A  o2  A  o2  (  o2  Z1  o2)  +  o2  Z'1  o2 


1 


detZ 


[1  -  (A  o2)2 ] Z: 


(3-26) 
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where 


(A  o2)2  =  A  a2A  o2 


^11^22  ^12  0 

0  ^11^22  ~  ^21 


(3-27) 


It  is  observed  that  Eq.  (3-27)  is  greatly  simplified  if  A  is  symmetric.  Therefore,  in 
this  thesis,  we  consider  only  the  case  when  A  =  A7  .  Then  A 12  =  A2,  so  that: 

(Ao2)2  =  (detA)/  (3-28) 


From  Eq.  (3-26), 


1  -  detA  N 
detZ  j 


(3-29) 


Eq.  (3-29)  can  be  satisfied  only  if 
Case  I 


1  -  det  A  ^ 
detZ  j 


=  ±  1 .  These  two  cases  are: 


Z  =  -  ZT  = 


0 


*12  *  0 


(3-30) 


detA  =  1  +  detZ  =  1  +  z12 


(3-31) 


Case  II 


Z  =  ZT 


(3-32) 
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detZ  +  detA  =  1 


(3-33) 


On  the  other  hand,  substituting  Eq.  (3-24)  into  Eq.  (3-22)  yields: 


Zll(A12_  A21^ 

(z12-z21)  An 

(3-34) 

^22^A12  ~  A21^ 

(zi2_z21)A22 

(3-35) 

A22  +  Z22A11 

2  Z21A12  =  ^12A21 

.(3-36) 

For  Case  1,  Eqs  (3-34)  and  (3-35)  require  that  An  =  A22  =  0 ,  Eq.  (3-36)  requires 
that  A12  =  A21  =  0.  Therefore  A  =  0.  From  Eq  (3-31),  1  +  z12  =  0.  Therefore  zn  =  +  i 
so  that  Z+=Z'=Z  =  ±o2. 

For  Case  II,  Eqs.  (3-34)  and  (3-35)  are  trivially  satisfied.  Eq  (3-36)  becomes: 

ZnA22  +  *22 An  -  2zi2A12  =  °  (3.37) 

Eq.  (3-33)  is  explicitly: 

*11^22  “  4  +  A11A22  “  A?2  =  1  -(3-38) 

The  sum  of  Eq.  (3-37)  with  Eq.(3-38)  is: 

( zu  +  An)  (z22  +  A22)  -  (z12  +  A12)2  =  det(Z  +  A)  =  detZ+  =  1  (3-39) 
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Subtracting  Eq.  (3-37)  from  Eq.  (3-38): 


(Zjj  ^22)  (^12  Aj2)  —  tiet(Z  A)  —  detZ  —  1  (3-40) 


In  summary,  two  sufficient  conditions  to  satisfy  Eqs.  (3-22)  and  (3-23)  have  been 
deduced  under  the  assumptions  that  A  is  symmetric  and  Z  is  invertible.  The  first  one  is: 

Z+  =  Z'  =  ±  o2  (3-41) 


The  other,  with  both  Z+  and  Z'  symmetric,  is: 

detZ+  =  detZ"  =  1  (3-42) 

It  should  be  noted  that  if  S  is  a  closed  surface,  the  impedance  boundary  condition 
closes  off  the  interior  of  the  surface  from  its  exterior.  Therefore  for  the  exterior  problem, 
det  Z +  =  1  is  sufficient  to  eliminate  the  on-axis  backscattering  from  a  body  of  90°  rotational 
symmetry.  Furthermore,  Z+  may  vary  with  location.  This  is  an  extension  of  Weston's  theory 
of  isotropic  absorbers  [4], 
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IV.  SCATTERING  OF  AN  ANISOTROPICALLY  COATED  CYLINDER 


In  this  chapter,  the  sum-difference  surface  current  formulation  developed  in  Chapter 
II  is  applied  to  the  problem  of  scattering  of  an  anisotropically  coated  tubular  circular  cylinder 
of  finite  length  and  negligible  wall  thickness.  Due  to  the  rotational  symmetry,  a  Fourier  series 
expansion  can  be  utilized  to  reduce  the  variables  of  the  problem  to  only  the  one  along  the 
axis  of  symmetry,  chosen  as  the  z-axis.  The  Fourier  components  Mn,  Nn  of  the  operators  M 

and  N  are  deduced  in  terms  of  the  Fourier  components  of  G(r-r0)  and  its  partial  derivatives. 
To  solve  the  set  of  integrodifferential  equations  so  obtained,  the  equations  and  the  surface 
currents  are  both  weighted  and  expanded  over  the  Chebyshev  polynomials.  The  resultant 
infinite  system  of  linear  equations  are  then  truncated  and  inverted  numerically. 

A.  GEOMETRY  AND  COORDINATE  SYSTEM 

Figure  4-1  shows  the  geometry  and  coordinate  system  of  a  tubular  circular  cylinder 
of  radius  a  and  length  2 h.  The  cylindrical  coordinate  system  (p,  <J),  z)  is  scaled  so  that  the 
surface  of  the  cylinder  S  is  specified  by  p  =  1  and  -  1  <  z  <  1 .  The  radial  vector  is  thus  given 
by: 


r  =  a  p  p  +  hzz 
=  a  p  (costjix  +  sin<j>y)  +  hzz 


(4-1) 


To  facilitate  representing  the  tangent  vectors  over  S  in  matrix  form,  we  chose  u  =  <J>  and 
v  =  z  as  the  orthonormal  basis  vectors  on  S,  so  that  n  =  p  is  the  outward  normal  to  S. 

To  properly  classify  the  polarization  of  the  incident  wave,  the  rectangular  coordinate 
system  (x,y,z)  is  determined  as  follows:  When  a  plane  wave  is  not  incident  along  the  z-axis, 
the  coordinates  are  chosen  so  that  the  wave  is  propagating  in  the  zx-plane  incident  from  the 
half-plane  containing  the  positive  x-axis.  Axial  incidences  then  are  considered  as  the  limiting 
cases  as  the  direction  of  incidence  approaches  the  positive  or  the  negative  z-axis  from  this 
half-plane.  Therefore,  even  for  axial  incidence,  a  linearly  polarized  incident  wave  having  its 
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Figure  4-1.  The  Geometry  and  Coordinate  System. 


electric  field  intensity  vector  E  pointing  in  the  y-direction  is  considered  a  TE  wave  while  one 

having  its  E  vector  in  the  zx-plane  is  considered  a  TM  wave.  In  the  spherical  coordinate 

system  (r,  0,<j)),  the  direction  of  propagation  of  the  incident  wave  k  is  given  by 

k  =  -  £  cos  0  .  -  x  sin  0  .  where  the  incident  angle  0  .  varies  over  the  range  0  <  0f  <  tz  .  For 

an  incident  plane  wave  of  unit  strength,  the  fields  are: 

TE  -  polariation: 

Emc  =  y  ei£mf 

Hinc  =  k  x  Einc  =  (f cos 6(.  -  zsinG,.)  eli'7 
TM  -  polarization: 
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Hmc  =  y  eli'7 

Einc  =  Hinc  x  k  =  -(icos 0(.  -  zsin 0;)  ei£'7 

where  k  =  kk.  Note  that  k*  r  =  l^cosB^  l2 p  sinSxostJ)  where  /,  =  kh,  and  l2  =  ka. 
On  the  surfaces  S,  the  tangential  components  of  TE-polarized  plane  wave  are: 
E%c  =  costj) 

=  -  cos0;sin<|)  el^‘r  (4-1) 

H‘nc  =  -  sin0;.e'*,f 


and  the  tangential  components  of  TM-polarized  plane  wave  are: 

E™  =  cosO^intj)  el*'r 

E™  =  sin  eig,‘*v  (4-2) 

H£c  =  coscj>e‘*,? 

B.  SCATTERED  FIELDS 

From  Eq.  (2-6),  the  scattered  fields  from  the  cylinder  are  given  in  terms  of  the  sum 
equivalent  electric  current  K  and  the  sum  magnetic  current  L : 

—E“(r)  =  -  -V  x  f'  f2’  L(r)G(r-r)d$dz 
k  I- Jo 

+  ‘  IX’  W-W.  (4-3) 

-  X IX' 


and 
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77 ■£“(*■>  =  7Vx  f'T  *<r«>  G{?-70)dbodzo 
^  ”1*'  0 

+  i  rr  ^ro]  G(-f~fo)d<&odzo 

J  ~\J0 

-.Lsj  nr2”  L(r  )  •  VG(r-f)4^ 
]^1  J-lJO 


(4-4) 


Their  components  in  cylindrical  coordinates  are: 
471 

Vz 


471  £pC(r)  =  7  /_[  f**  ( $  *  4>0)  ^  (ra)  G  (r-ro)  d<\>0  dz0 


i  ^ 

-  — - — —  f1  f2n  L  (r  )  G(r-r  )  d( J) 

4p  94>  J- 1  Jo  z  0  000 

+  f  h  C  (p‘^}  K*(V  G(p~°  dz<> 

+ 1  /-,  r 

+  7^afk/-r  fWOF-WA. 


(4-5) 


4  7t 


=  | /;  fwwop-v*.*. 


/,  dz 
1  5  /•  1  f2it 


+ - f 1  r2n  L  (r)  G(r-r  )  4(|>  4z 

/2  5p  J-i  Jo  z  “  °  0  0 

+  '  /.j  .C*  ^  G(F“ Fo)  <*4>„ 

*pf,rw^ 

*i£tmf-‘CJC‘(oa(’’'°d*'‘K 


(4-6) 
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E*c(r)  =  — - — —  f 1  [2*  ($•  <j)  )EA(r  )  G  (r-r  )  d§  dz 
Z,/2  *  l2 p  0p  J- 1  Jo  +V  °  V  o'  Vo  o 

+  — — —  f 1  f2n(P  •$  )  L.  (r  )  G(r-r  )d'(j)  <iz 
Z2P  0<J)  J-l  Jo  °  *  0  o’  Va  0 


imfXK*(F’)0(F-vd*°dz° 


i  1+  —  —  f  f  K7(r)G(r-r)d<t>dz0 
f  dz2  J-l  Jo  z  °  o  0  ° 


—  =  —  f1  f2”($*<i>  )**(?,)  G(r-rJ<*b  & 

/j/2  p  Z,  dzJ-iJo  °  *  0  °  0  ° 


+  — — —  f 1  f2K  K  (r  )  G(f-f )  dz 

Z2p  0cj)  J-i  Jo  z  0  0  0  0 

+ '  f-i  f02K  ( ^  *  ^o)  L* (Fo)  G(?"  ^  *0 

+  —  — —  fl  f2nLs(r„)G(r-r)d<t>  dz 

;2  0p0(j)  J-l  Jo  *  0  0  0  0 


i  d 2  ri  f2n 


/,Z2  3pdz  J -l  Jo 


L(nG(F-F)4^ 


t*;C(rl=  i  I/-!  fv-tww-v*.*. 

-  t4z  f!  f.2* 


Z2  0p  J-i  Jo  z'  °'  '  °  0 

+  i  [l  C\$-b0)L^r0)G(?-ro)d$0dz0 

J  -1  J  o 
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4tT  sc , 


-kf', 


hh 


1 i&Cfww-v*.*. 


(4-10) 


+ i 


1  + 


1JL 

if  dz2 


/  C*  ^)G(F-F)  d$odz„ 

J  - 1  j  o 


Note  that  in  the  above  equations, 

P  =  p0cos((j)-({)rt)  +  4>0sin((j)-(f)o)  (4-11) 

I  =  _  P0sin((J)  -c])^)  +  4>0cos  (4>  -  cj>0)  (4-12) 


Because  of  the  rotational  symmetry  G(r-r0 )  depends  on  4>  -  4>o  and  Fourier  series 

can  be  introduced  to  eliminate  the  variable  (J) .  Define  the  Fourier  series  expansion  of  a 
function/(4>,z)  by: 

00 

/(<M)  =  £  ein*fn(z)  (4-13) 

n  =  -oo 


then 


ik\r-rQ\  oo 

=  E  e^GMz-z,  | ,  4,  p) 


k\r-fo\ 


(4-14) 


and 


30 


(4-15) 


J  -71  ZTC 


Eqs.  (4-5)  to  (4-10)  become: 

2  j — i  sc  s  _  ^  1  d 


-  %  W, 

-  t  /.:  wa 

4  W  I<v, W 

+  _L  Af 

/j/2  0pJ-l 


^(p’z)  ‘  4  I/.',  W<v.-<v,)^  /.;  we. A 


+  i 


/.;  w 

n  />i  |  5 
Z,/2p  J-i 


rr(G.-i+G.tI)-i-G. 

2  /2P 


G  Jz 


4£"<p,z)  =  2ip  /-', 


/',  w,^ia 


/»  1 

/- 

1 

kw] 

g,a +i  L[  w 


(4-16) 


(4-17) 


(4-18) 
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1,1  H‘"'  P,Z>  21,  dzf- 1 


1  n 


f_[  w  lG..1-GH.i]dz0-±-±f_,i  VO 

Z-,  * 

i 

Z]/2  dpJ-i 


(4-19) 


G„<fe„ 


*Op.z>  =  -  jrfj',  ~  f[  WG,dz, 

*•  /'  VO 

J’1  Z  /oP 


v2 


a?  /.;  t|wi 


(4-20) 


V2 


«P.d  --  /'  A*,^[(n-l)C..1-(»*l)G..i]<fe. 

is/>eA 
4|  /_;  t^-wi  ga  ♦<  /.;  w 


(4-21) 


Note  that  in  Eqs.  (4-16)  to  (4-21),  Gn  stands  for  Gn(Z,  |z-zj  ,/2,  p) .  In  shifting  the  z-derivative 

to  the  current  densities  A^(z0)  and  LJ^Zq)  ,  the  fact  that  edge  conditions  [5]  require  that  these 
components  of  the  scattering  currents  vanish  as  |  z0|  approach  1  from  below  is  used. 
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As  p  -  1*,  the  operators  Mn  and  Nn  can  be  deduced  from  Eqs.  (4-16)  to  (4-21): 


(4-22) 


N‘-"sikTz 


(4-23) 


Note  that: 

•  (  - 


U 


\ 


\  dp  p  =  r  +  aPlp=r, 


\\z  20 1 ,  l2,  p )  ^  Gn  ( /j  |  z  Zq  | >  l2, 1 )  (4-24) 


Assuming  that  Z  is  invertible,  from  Eq.  (2-23),  the  equations  for  the  sum  currents  are 
therefore:  - 


M 

n 

-N  1 

n 

Kn 

Kn 

Emc 

tan,  n 

+  R 

=  2 

Nn 

K 

Ln 

L 

H‘nc 

n 

tan,  n 

(4-25) 


where 
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(4-26) 


R  = 


Z-AZ1  A 
-z'a2Z-1  A 


-i  AZ  1  a2 
°2Z~lo2 


and  the  vectors^  ,  L  ,  E,mc  ,  L,mc  are  two  dimensional  column  matrix  representations  of 

n’  tan,n’  tan,n  A 

the  respective  tangential  vector  fields  on  S  over  the  orthonomal  basis  <f>,  z-  The  incident 
fields  are,  for  TE  and  TM  polarizations  respectively: 

TE: 


C  =  (-  0”’  V«(4  sin  6j) e  ~il}Zcos6‘ 


H, 


n(-i)n  cos0. 


4>,w 


l2  sin  0; 


/n(/2sin0(.)e 


-/7Izcos0i 


n-\j  n  \  »  -i^zcosOj 


tfz,„  =  -sin0(.(-i)n~Vn(/2sin0(.)e 


(4-27) 


TM: 


£"!  =  J'(l2 sine,) e-“'zm6i 


/2sin0(. 


E™  =  sin  0(. ( - i)n  Jn ( l2  sin  0;.)  e  ~il'lcosQ‘ 
Kn  =  (-0n-lj'(l2  sin0.)c"a,zcos6i 


(4-28) 


C.  TRANSFORM  TO  SYSTEM  OF  LINEAR  EQUATIONS 

Since  the  surface  current  components  K^,z0)  =  0(l -Zg)’1-  and 

Kz($,z0)=0(  1  -z02)1/2  as  |zj  -  1“,  representations  of  K^n(zo)  and in  Chebyshev 

aZo 

polynomials  of  the  first  kind  combined  with  the  weighting  factor  conform  to  the  proper  edge 
behavior  of  the  currents: 
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(4-29) 


K*,„(z„)  -  —  E  Kl„  T  (Za) 

TUSinV  q=0  ^  1 

=  - —  E  Kln  COS<?V 

ttsinv  o 

K„fe„)  -  1  E  sin(<?+l)v 
^  ^  “0 


where  rp(zo)is  the  Chebyshev  polynomials  of  the  first  kind  and  z0  =  cosv,  - 1  <  zQ  <  1 . 
Similarly, 


L^,SZo}  =  — —  J2Li,n  C0S4V 

T  TTSinv  <7=0 

Lz,n(zo)  3  —  E  Lz,n  sin(^  +  l)v 

^  4=0 


(4-30) 


For  an  invertible  Z,  the  coefficients  in  the  above  equations  are  to  be  determined  from 
Eq.  (4-23).  To  make  use  of  the  orthogonal  property  of  the  Chebyshev  polynomials,  the  factor 
sin  v  sin  (p  +  1)  vfor  p  >  0  is  multiplied  to  both  sides  of  Eq.  (4-23)  before  an  integration  over 
the  range  of  v  is  carried  out.  The  results  are  described  term  by  term  in  the  subsection  to 
follow.  This  procedure  creates  an  infinite  system  of  linear  equations  to  be  solved  numerically 
after  it  is  truncated  at  an  appropriate  order  determined  by  the  electrical  size  of  the  cylinder. 


1.  Incident  Fields 


ffinc,p 

tan,n 

2 

rn  dv 

sinv  sin(p  +  l)v 

Einc 

tan.n 

frinc,  p 

P  +  1* 

1 o  tu 

Hinc 

tan,n 

tan,n 

(4-31) 
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TE: 


Ei;cf  =  (-on+P-1(7p_1(/1cose,.)+7p+1(/1cose,.))7,;(/2sine,.) 

=  -  y/'ife"  '/p+i(Zicos9')  Jr»(/2sin0‘)  (4_32) 

-  -(-i)n+psinei.(7/,.iaicosei.)+7p+1(/1cose,.))7n(/2sin0l.) 

TM: 

-  2.W.(~.^  ^+,(^OS0(.)  7n(/2sin0,) 

33) 

Einc,P  =  (_On^sin0[  (  ^^cosO,.)  +  7p + , (/, cos 0,) )  7n(/2sin  0;.) 

<c/  -  (-0n+p_1  (  /p.j^cosO,.)  +  Jp + ! (/, cos 0;) )  7^2sin0;) 


where  7 '( • )  is  a  derivative  with  respect  to  the  argument.  Note  that  as  0 .  approaches  0  or  7t , 
only  n  =  ±1  terms  are  nonzero.  Hence  only  n  =  ±1  currents  exist.  For  axial  incident  when 
0,  =  0: 

TE: 


77  inc,p 

^.l 


inc,p 


_  (“*)P  r  r/N  =  17 inc,p 
- : —  7  -  ^,-1 


nr 


inc,p 


K,  i(V  = 


(4-34) 


TM: 


inc,p  (O^^t  /  /  \  _  r  'nc'P 

i<M  =  - -  7p+1(/i)  -  -%-i 


Yjinc,p  _ 


(-O' 


Jp-iC'i) 


Trine, p 


(4-35) 
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2. 


The  R  -  Matrix  Term 


—^r  fU—  sinv  sin(p + 1  )vK.  (cos  v) 
p  +  lJO  K  V 

=  - Y  Kin  rdv  sin(p  +  l)v  cosgv 

( p+l)n 2  <7=o  Jo 

- 1  K"kI, 

9=  0 


where 


4P’9  - 
a4>  " 


0 

4 


[  n2(p+q  +  l)  (p-q  +  l) 


p+q  odd 
p+q  even 


_2_r*dv  s-nv  sjn(p+j)v  g  (cosv) 

p+lJo  71 

2 

=  -  Y  K?  fndv  sinv  sin(p+l)v  sin(#+l)v 

(p+ 1)tt2  q= o  ’  Jo 

-  ixMJ5i 

,.o 


where 


0 

_ -8(g+ 1) _ 

n:2(p+^  +  l)  (p-g+l)  (p+q +3)  (p-q- 1) 


/?+<?  odd 
/?+g  even 


Therefore: 


(4-36) 


(4-37) 


(4-38) 


(4-39) 
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p+ 


2  nn  dv  .  .  ,  t .  n  I  n 

-/  —  sinv  sin(p+l)v  R 
1  Jo  ft 


K 


■  E  * 

<?  =  0 


=  E 

9=o 


Ar‘ ¥  ooo 


0  A['q  0  0 


0  0  AZ'q  0 


0  0  0  A['q 


K‘ 


L 


K 


(4-40) 


where  Rp,q  is  a  four-by-four  matrix,  given  in  terms  of  R  by 


Rfiq  =  Ru  Kq 

R£'  »  Ra  A™ 

Rf?  -  R„  K" 

R,?  =  Ru  A™ 


f  for  i  =  1,  2,  3,  4 


(4-41) 


Note  that  Rf}'q  =  0  if  p+q  is  odd. 

3.  The  Mn,  Nn  Operators 
Define  the  4x4  matrix X„‘ q  by: 


-fn—  sinv  sin(p  +  l)v 
1  Jo  ft 


\M  -Nn] 

Kr 

oo 

n  n 

n 

‘  E  X ,M 

n 

N  Mn 

n  n 

»-o 

T  c 

(4-42) 


2  /o’v  Sinv  Sin(p*1)v  M »■" 


P+1 


.  dv„ 


=  r^[cospvcos(p+2)v]r^cosqvU(Gn_1+Gn+1)-^-Gn 

q=Q  /?+!  *^0  TC  JO  U  2  ^ 

oo 

EVA9  ]7? 

An,ll  A(J ),n 

9  “0 


K. 
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and, 


YP’<1  . 
An,  11  " 


i/j/2 

T* T 


^  (ap'q  ~ap+2'q  +  rzp'q  np+2'q\-  n  (ap’q  -ap+2'q 

T^n-l  °Vl  +  0r/i  +  l  °Vi+l  ' 

2  /Z 

to 


(4-43) 


^-/;f  sinv  sin(p+l)v  Mn  21  [^,n(z0)] 
-  2n  rndv  . 


<iv 


E- Ln  rnciv  .  ,  a  ^  ^  ^ 

— -  —  sin(p  +  l)v—  - cosgv  G 

q=o  p+ 1  Jo  it  dvJo  7t 

=  i,2n  [*—COS(P  +  1)v[%—COSqVoGnK£n 

q=0  J0  71  Jo  n 

-  r  ? 

“  An,21  A4),/i 

<?=0 


and, 


X«  =  2nGr'-‘' 


(4-44) 


2  rndv 


p+ i 


f ~  sinv  sin(p  +  l)v  7V„  u  n(z0)] 

^0  71 


*=0  P  +  1  ^0  71  8W0  K 


<iv 


oo 

V?'?  jy9 
Z^An,  31  A<t >,n 
9  =0 


and, 


(4-45) 
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~ -/;f  sinv  sin(p  +  l)v  N  21  [^n(z0)] 

=  — —  f*—  sinvsin(p+l)v  fn — -cosqv 

q=o  p+ 1  Jo  n  J  on  0 

(»-  l)GJI_1  -(n  +  l)Gn+. 


2  a/. 


(G 


„-l+G*  +  l) 


F? 

a4>,« 


-  V  yP'%  y^ 

A Lu  An,  41  A(J>, 
4=0 


and, 


Yp’q 

An,41 


2(p+l) 


(n  -  iKG"-G,f:l2',i-(n  *  ikg;;’  -g,'’;,2'') 

- |  J-  cor,  -  g::,2-’ + g';’  -  G„'’,f<) 


(4-46) 


and. 


^t/o’tT  sinv  sin(,,*1)v  M“'12  [^”a”>1 

■  X  [’— sinvsin(p+l)vP — -(«  +  l)cos(g  +  l)v  5^' 

9  =  0  p  +  1  Jo  1  Jo  71 

CO 

-  yP'Q  Y? 

2-j  A«,  12  Az,rc 
?=0 


yP^ 

a/i,12 


W(£+1)  r/^P>4  +  1  ^fP+2,^  +  1-. 

”™  -  L  -1 

p+i 


(4-47) 
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sinv  sin(p+1)v  K22 


.dv 


Y - -  I  f71—  sinvsin(p  +  l)v  fn — -sinv  sin(<?  +  l)v  G 

q =0  P  +  1  [JO  TC  JO  71 

+  J_  f  sin(p  +  l)v  JL  (4  +  l)cos(g  +  l)v  G  l  K.q 

1 2  Jo  7C  dvJo  n 

L\ 


S'  yM  v q 

Z-j  A/z,22  A2,h 

q= 0 


and, 


Yp'q  - 
A«,22  ~ 


i7j  /2 


+  QP  +  2W  +2  _ +2  _  QP  +  2^\  _  4(£^2)  qP  +  M  + 1 

.  1  '  w  rc  /x  '  o  n 


p+i 


h 


'(4-48) 


^r/o’^  Si"V  sin^"'>v^,2  [^.(z,)] 


.  dv 


=  5^-^-  f71— sinvsin(p  +  l)v  fn — -  sinv  sin  (q  +  l)v 
q= 0  p+ 1  Jo  7T  Jo  it 


3^-1^ 


EvP.?  js<l 
An,32  Kz,n 

q=0 


and, 


Yp'q  - 
Anf  32  " 


/i/2 


4(p  +  l) 


—  [G*« 
di2  n 


QP*2,q*2  _Qp*l,q  _Qp,q *2- 


Note  that,  for  n=0,  the  expressions  simplify  to: 


(4-49) 
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il,  L 


Ku  =  — 4-[Gf’9-Gf+2’9] 

P  + 1 


y/>»<7  _ 

Ao,41  “ 


vP>?  _ 

A<?,22  “ 


YP'V  - 
Ao,  32 


Z|  [  (Gf  ’ 9  -  Gf +2'9)  +  i—  (Gf’9  -  Gf + 2'9)] 


P  +  1 
/  /]  /2 

~T” 

hh 


2  5L 


1 _ +  _  4(g  +  l)^.p  +  i,iy  +  i 

t  v  '~7o  'Jo  '-'o  ^0  /  o 


P  +  1 


K 


l',  Q 


(4-50) 


4(p  + 1)  dl 


(qP'Q  +  qP +2'Q +2  _  qP’I  *2  _  qp  *2-^ 


yP’Q  -  Y p,q  -  Yp'q  —  Yp,q  —  A 
Ao,  21  _  ao,31  _  ao,44  -  Ao,  12  _  u 


Also  note  that  Xf’9=  0  if  p+q  is  odd. 


From  the  symmetry  of 


M 

-N 

n 

n 

K 

we  can  deduce  that: 


YP'Q  - 
Atz,  13  " 

_  YP'q 

Ati,31 

YP>q  - 
Any  14  " 

An,  32 

YP'1  - 
Any  23  “ 

YP'q 

A71, 41 

YP'q  - 

An,  24  ~ 

y/?>4  _ 

Ati,42  “ 

YP'Q  - 
An,33  ~ 

YP'q 

An,  11 

VP*?  - 
ati,34  " 

Yp'q 

An,  12 

Yp'q  - 
A /i,43  “ 

YP'q 

Anyl\ 

yP.fl  _ 
A 71,44  “ 

YP'q 

Any  22 

(4-51) 


Note  that  for  p+q  odd, 


vP.9  _  Yp-q  -  Yp'q  -  Yp'q  -  Yp'q  -  Yp  q  -  Yp’q  -  Yp'q  -  A 
An,  1 1  _  An,41  “  An,22  “  An,32  "  An,23  “  An,33  "  An,14  "  An,44  "  U 


(4-52) 


and  for  p+q  even, 


Yp’q 

An,21 


Yp’q 

An,31 


yP>9  _  yP'Q 
An,  12  "  An,  13 


_  yP'Q  -  v/7’^  -  n 

“  Atz,43  “  Ati, 34  “  U 


(4-53) 


The  integrodifferential  Eq.  (4-25)  is  thus  transformed  into  the  infinite  system  of  linear 
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equations  of  the  unknown  coefficients 


[Xpn'q  +  R™] 


ginc.p 

n 

=  2 

tan, /i 

V 

ffinc.p 

n  j 

tan,n 

(4-54) 


which  is  to  be  solved  numerically. 

D.  RADIATION  IN  THE  FAR  FIELD 

In  the  far  field,  we  can  write  Eq.  (4-3)  as: 

77- E  sc(r)  ~+i  [idzo[ 2nd$o [ L(rg) xr)G(f-ro) 

+  i  f  ldzof  27I#0  {  ^  (ro)  -  if  [  K^(ro)  sin  0  sin(4>  -  <{)„) 

+  Kz(ro)cosQ]}G(r-ro) 


(4-55) 


or  equivalently, 
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=  ~y  j'dz0  /;  {  ©  [  ^  cos  0  sin(cj)  -<$>o)  -Kz  sin  0  +L^cos(<$>  -<$>o)] 

+  [  A^cos^  -<{)o)  -L^cosG  sin(cj)  -4>o)  +Lzsin0]  |  e  ‘4Mocose+/2sm0cos(|l>-<tv 

=  !4r  E  (-ov*+  £  (-o'- 

Z  fl  =  -00  p  =0 


L  V 


i7Asi„e)^®  -y'(/2sin0) Li„ 

z 


y  (Zjcos0) 


(4-56) 


7n(/2sin  0) sin  6  Zi cos  0)  +  ^+2(zi  cos  0)) 


-  <$» 


7'(/2sin  6)  +  iV„(/2sin  0)  L 


7  (Zj  cos  0) 


(Z2  sin  0)  sin  0  LzPn  ( 7p  (/,  cos  0)  +  7p+2(Z,cos0) ) 


As  0  -^0 ,  only  the  n  =  ±  1  terms  are  nonzero  in  this  limit,  and 
p  =  x  coscj)  +ysin<|),  <$>  =  —  jc  sin  cj>  +  ycos<|).  We  have  &  ±  z<|)  =  (jc  ±  y)eri<^: 


Esc(r)_  V 2 
G(f)  "  4 


E  ( -  o '  { *  [  (*£  i  -  K  - , ) + [Ll . +  ^  - 1 )] 

p=0 


(4-57) 


Similarly,  as  0  --  it  ,  only  the  n  =  ±  1  terms  are  nonzero.  0  =  -  p  =  -icosc{)  -  ysincj), 
$  =  -xsin(J)  +  y  cos (J> .  We  have  @  ±  z<|)  =  -(x  +  y)e±,<t>; 


T%  "Wte -k*.-‘)-‘K  *Li-\ 


(4-58) 
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E.  INSIDE  AND  OUTSIDE  SURFACE  CURRENTS 


From  Eq.  (2-20)  we  have: 

KP~KP  If  Z“‘A  iZ~lo2  1  \kp 

L+p  -  Lv  -io2(Z-AZ-‘A)  -o,AZ-lo, 

n  n  L  z  z  zj  /i 


(4-59) 


Therefore  the  following  matrix  equation  is  obtained: 

Cl  J  f  ^A  iZ-‘o2  Ilk 
CJ "  2I  +  h'^cz-Az-^)  -o2az_i o2j  J  [r; 


(4-60) 
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V.  COMPUTATION  AND  RESULTS 


The  solution  to  the  problem  of  the  scattering  of  an  anisotropically  coated  tubular 
cylinder  of  finite  length  as  formulated  in  the  previous  chapter  has  been  coded  in  FORTRAN 
and  tested.  The  program  listings  are  included  in  the  Appendix.  Computation  has  been  carried 
out  on  the  32-bit  Sun  SPARC  Station  running  under  the  Unix  operating  system  in  the 
Electrical  and  Computer  Engineering  Department  and  the  Computer  Center.  The  evaluation 
of  the  double  series  expansion  coefficients  of  the  Green's  function  and  its  derivatives  for 
greater  values  of  ka  and  kh  have  also  been  done  on  the  64-bit  Cray  Y-MP  EL98  in  the 
Visualization  Lab  so  that  the  accuracy  of  the  results  can  be  accessed. 

The  program  accepts  the  geometry  of  the  cylinder,  the  surface  impedance  matrices, 
the  incident  wave  frequency,  direction  and  its  polarization  in  terms  of  TE,  TM  or  some  linear 
combination  of  the  two.  It  computes  the  sum  surface  currents  and  the  far  field  radiation  in 
the  directions  specified,  including  the  strength  and  the  phase  of  each  field  component.  It  also 
breaks  down  the  sum  surface  currents  into  the  outside  and  the  inside  currents. 

In  this  chapter,  some  interesting  results  of  computation  for  the  scattering  of  a  tubular 
cylinder  having  the  length-to-diameter  ratio  h/a  of  4  or  6  are  presented.  All  figures  are 
attached  at  the  end  of  this  chapter.  The  wave  is  incident  from  the  positive  z-axis  and  is 
polarized  in  the  y  direction. 

A.  COMPARISON  WITH  EXPERIMENTAL  DATA 

For  backscattering  from  a  perfectly  conducting  tubular  cylinder  of  a  y-polarized  plane 
wave  incident  from  the  positive  z-axis,  experimental  data  are  available  [6]  over  a  frequency 
band  of  well  beyond  three  octaves,  with  the  circumference-to-wavelength  ratio 
ka  =  2rt<3/ X  varied  from  0.9448  to  3.3152.  These  data  are  measured  with  two  sets  of  four 
cylinders  each;  one  set  having  h/a  =  4,  the  other  with  h/a  =  6.  Both  data  sets  use  the  inner 
radii  of  the  tubular  cylinders  as  the  parameter  a  [7].  The  cylinder  length-to- wavelength  ratio 
2 h!  X  varies  from  1 .2030  to  4.2210  for  the  h/a  =  4  case  and  from  1.8045  to  6.33 15  for  the  h/a 
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=  6  case. 

The  experimental  data  are  plotted  against  the  results  of  theoretical  computation. 
Figure  5-1  shows  the  backscattering  from  the  set  of  cylinders  having  h/a  =  4.  Figure  5-2 
shows  those  data  from  the  set  of  cylinders  with  h/a  =  6.  The  output  of  theoretical 
computation  follows  the  measured  data  very  closely.  Note  that  the  cutoff  frequency  of  the 
dominant  circular  waveguide  mode,  TE,„  occurs  at  2  hi  X  =  2.344  for  the  cylinders  with  hi  a 
=  4  and  at  2 h/X  =  3.516  for  those  with  h/a  =  6. 

B.  NULL  ON-AXIS  BACKSCATTERING 

Three  different  ways  of  coating  a  surface  impedance  Zs  on  a  tubular  cylinder  of  h/a 
=  6  are  considered:  Both  on  the  outside  surface  and  the  inside  surface;  Only  on  the  inside  and 
leaving  the  outside  surface  perfectly  conducting;  Only  on  the  outside  and  leaving  the  inside 
surface  perfectly  conducting.  Here  Zs  has  the  elements  zsU  =  0.5,  zs22  varies  from  0.1  to  5,  zsl2 
=  zs2i  =  0.  At  a  fixed  frequency  for  which  2h/X  =3.194,  slightly  below  the  T£  circular 
waveguide  dominant  mode  cutoff  of  3.516,  the  scattered  fields  are  plotted  in  Figures  5-3 
through  5-5.  Figure  5-3  shows  the  results  of  computation  for  the  case  Z+  =  Z  =  Zs  =  Z.  As 
22  is  varied  through  2,  the  backscattering  cross  section  vanishes  as  predicted  in  Chapter  3. 
Figure  5-4  shows  the  results  for  the  case  Z+  =  0  and  T  -  Zs.  As  the  impedance  on  the  inside 
surface  is  increased,  the  excited  field  inside  the  tubular  cylinder  is  dissipated  and  the 
backscattered  power  decreases  exponentially.  Figure  5-5  shows  the  results  for  the  case  Z+  = 
Zs  and  Z  =  0.  The  backscattered  power  drops  off  rapidly  at  first  as  the  impedance  on  the 
outside  surface  is  increased.  But  the  cross  section  quickly  settles  down  to  a  fixed  value 
presumably  due  to  the  current  excited  on  the  perfectly  conducting  inside  surface  of  the 
cylinder. 

Results  of  computation  for  the  same  configurations  but  at  a  higher  frequency  for 
which  2 h/X  =  4.865 ,  above  the  TEn  circular  waveguide  dominant  mode  cutoff  of  3.516,  are 
plotted  in  Figures  5-6  through  5-8.  Now  that  the  incident  wave  can  propagate  through  the 
cylinder  in  the  dominant  waveguide  mode,  the  backscattering  cross  sections  are  about  an 
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order  of  magnitude  smaller  than  in  previous  cases.  Figure  5-6  shows  the  results  of 
computation  for  the  case  Z+  =  Z"  =  Zs  =  Z.  Again  the  backscattering  cross  section  vanishes 
as  z22  is  varied  through  2.  Figure  5-7  shows  the  results  for  the  case  Z+  =  0  and  Z'  =  Zs  while 
Figure  5-8  shows  the  results  for  the  case  Z+  =  Zs  and  Z‘  =  0.  It  appears  that,  above  cutoff,  the 
contribution  to  the  backscattering  cross  section  from  the  inside  of  the  tubular  cylinder  is 
minimal:  once  the  outside  current  is  reduced  by  the  increase  in  surface  impedance,  the 
backscattering  cross  section  is  reduce  by  more  than  10  dB  as  shown  in  Figure  5-8.  When  the 
impedance  coating  is  applied  in  the  inside  surface,  the  maximal  reduction  in  the 
backscattering  cross  section  is  only  about  1.2  dB. 

C.  FREQUENCY  DEPENDENCE 

The  axial  backscattering  of  the  two  cases  when  the  tubular  cylinder  is  coated  only  on 
the  inside  or  only  on  the  outside  with  zs!1  =  0.5  and  %22  =  2  are  investigated  for  different 
frequencies  with  2h/X  varying  from  0.1  to  7.5.  Figure  5-9  shows  the  results  of  the  case  when 
only  the  inside  surface  is  coated  so  that  the  backscattering  is  mainly  due  to  the  current 
excited  on  the  outside  surface.  The  reflection  from  the  ends  of  the  cylinder  causes  the 
fluctuation  in  backscattering  cross  section.  Being  waves  in  free  space  on  the  outside  of  the 
cylinder,  the  maxima  and  minima  are  evenly  spaced  with  the  minima  occurring  when  2 hi  X 
is  a  multiple  of  half  integer.  Figure  5-10  shows  the  results  when  only  the  outside  surface  is 
coated  and  the  current  on  the  inside  surface  dominates  the  contribution.  The  distinct  feature 
in  this  case  is  that  the  backscattering  cross  section  does  not  fluctuate  with  varying  frequency 
below  the  waveguide  mode  cutoff.  The  incident  wave  is  able  to  penetrate  deeper  into  the 
cylinder  with  increasing  frequency,  resulting  in  a  constantly  rising  strength  of  the 
backscattered  field.  Once  beyond  the  circular  waveguide  mode  cutoff,  the  wave  can  pass 
through  the  cylinder  in  the  TEU  mode  and  the  backscattering  diminishes.  The  oscillation  in 
the  cross  section  at  these  higher  frequencies  represents  the  interference  of  reflected  waves 
at  the  ends  of  the  tube  and  the  separation  between  maxima  and  minima  should  be  determined 
by  the  guide  wavelength  at  the  particular  frequency.  These  two  situations  should  be 
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compared  to  Figure  5-11  which  shows  the  results  when  both  sides  of  the  cylinder  are 
perfectly  conducting  and  the  current  can  flow  freely.  The  distinct  notch  in  the  cross  section 
near  the  TEn  mode  cutoff  at  2 hi  X  =3.516  and  the  subsequent  faster  variation  in  the  cross 
section  shows  the  combination  of  the  two  distinct  features  of  Figures  5-9  and  5-10. 

D.  COMPUTATION  ACCURACY 

The  main  difficulty  encountered  in  the  computation  is  the  evaluation  of  G„q(lvl2) 
and  its  /2-derivative  by  double  power  series  sum  when  becomes  large.  Computations  for 
Figure  5-11, 5-13  and  5-15  use  the  G%'q  values  evaluated  with  the  Cray  computer  which  has 
a  128-bit  double  precision  number.  Computations  for  Figures  5-12,  5-14  and  5-16  use  the 
G%'q  values  evaluated  with  the  Sun  SPARC  Station  which  has  a  64-bit  double  precision 
number.  For  2 h/X  greater  than  about  6.2,  the  SPARC  Station  fails  to  provide  accurate 
results. 

Figures  5-13  and  5-15  are  axial  backscattering  from  a  cylinder  of  h/a  =  6  coated  with 
the  impedances  having  the  elements  zn+  =  z12  =  0.5,  z22+  =  Z\\  =  0.4,  zn+  =  z2l+  =  Zl2~  =  z2{ 
=  0.3.  Figure  5-13  shows  the  co-polarized  backscattered  field  while  Figure  5-15  shows  the 
cross-polarized  backscattering. 
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VI.  CONCLUSIONS 


In  this  thesis,  the  sum-difference  surface  current  formulation  is  introduced  for  solving 
electromagnetic  boundary  value  problems  when  impedances  are  specified  on  both  sides  of 
a  surface.  For  an  impedance  coated  body,  the  body  can  be  treated  as  being  a  surface 
separating  the  space  into  two  regions  of  identical  medium.  For  an  exterior  problem,  the 
impedance  normalized  to  the  medium  on  the  inside  surface,  Z  ~ ,  can  be  chosen  arbitrarily; 
and  for  an  interior  problem,  that  on  the  outside  surface,  Z + ,  can  be  arbitrary.  The  choice 
when  Z "  =  -  Z +  is  of  particular  interest  because  the  integrodifferential  equation  has  only  the 
sum  of  the  equivalent  electric  surface  currents  on  the  outside  and  the  inside  surfaces  as  its 
unknown  to  be  solved. 

This  formulation  preserves  the  duality  nature  of  Maxwell’s  equations  and  carries  it 
over  into  the  algebraic  form  of  the  integrodifferential  operators  in  the  equations  for  the  sum 
currents.  Since  a  90°  rotation  is  equivalent  to  undergoing  a  duality  transform  for  an  incident 
plane  wave,  this  particular  symmetry  in  the  algebraic  form  of  the  operators  leads  to  the 
sufficient  conditions  that  if  Z+  =  Z'  =  ±o2,  or  if  Z+  and  Z'  are  symmetric  and 
det  Z  *  =  det  Z  “  =  1  ,  the  on-axis  backscattering  of  an  anisotropic  impedance  coated  scatterer 
having  a  90°  rotational  symmetry  will  be  eliminated.  Note  that  in  the  symmetric  case,  Z + 
and  Z '  may  vary  with  location.  This  is  an  extension  of  Weston's  result  [4]  for  which  the 
surface  impedance  is  isotropic. 

A  FORTRAN  program  has  been  written  which  accepts  the  geometry  of  the  cylinder, 
the  surface  impedance  matrices,  the  incident  wave  frequency,  direction  and  its  polarization 
in  terms  of  TE,  TM  or  some  linear  combination  of  the  two.  It  computes  the  sum  surface 
currents  and  the  far  field  radiation  in  the  directions  specified,  including  the  strength  and  the 
phase  of  each  field  component.  It  also  breaks  down  the  sum  surface  currents  into  the  outside 
and  the  inside  currents. 

The  results  of  computation  using  this  program  agree  with  measured  data  of 
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backscattering  from  conducting  tubular  cylinders  over  a  frequency  band  of  more  than  three 
octaves.  For  a  cylinder  coated  with  surface  impedance  matrices  satisfying  the  criteria  for  null 
on-axis  backscattering,  the  numerical  computation  also  validated  the  theoretical  assertion. 

Difficulties  have  been  encountered  about  the  computational  accuracy  in  the 
evaluation  of  the  double  series  Chebyshev  expansion  coefficients  of  the  Green's  function 
G^'q(l{,l2)  and  its  /2-derivative  by  double  power  series  sum  when  the  length  of  the  cylinder 
is  large  compared  to  the  wavelength:  a  compiler  with  a  64-bit  double  precision  number  can 
only  handle  a  cylinder  having  a  length  up  to  about  6.2  wavelengths.  Further  work  to  explore 
the  feasibility  of  asymptotic  evaluation  of  these  coefficients  is  recommended. 
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APPENDIX  PROGRAM  LISTING 

A.  INCLUDE  FILES 


REALTP . INC 

C  TYPE  STATEMENTS  FOR  REAL  AND  INTEGERS  AND  DEFINITIONS  OF  CONSTANTS. 
IMPLICIT  DOUBLE  PRECISION  (A,  B,  D-H,  O-Z) 

IMPLICIT  INTEGER* 4  (I-N) 

PARAMETER  (PI=3 . 141592 653 58979323 8462 6433 8327D0 , PI2=PI+PI , 

+  PISQ=PI*PI) 

PARAMETER  ( ONE=l . DO , TWO=2 . DO , THR=3 . DO , HXD=1 6 . DO , ZERO=0 . DO  , 

+  DEGPI=180. DO, EPS8=2. 22044 6049250313D-16) 

PARAMETER  (ONEN=-ONE,  HALF=ONE/TWO ,  THIR=ONE/THR,  QUAR=HALF*HALF) 


CMPXPT . INC 

C  IMPLICIT  TYPE  STATEMENT  AND  CONSTANTS  FOR  COMPLEX  NUMBERS. 
IMPLICIT  DOUBLE  COMPLEX  (C) 

PARAMETER  (CZERO= ( 0 . dO , 0 . dO ) # CONE= {I . dO , 0 . dO) ) 

PARAMETER  (CONEN= ( -1 . dO , 0 . dO ) , CI1= { 0 . dO , 1 . dO ) , CI2= (0 . dO , -1 . dO )  ) 


B.  INPUT  DATA  FILES 


CLYGEOM .  PRM 

30  Maximum  kh  (integer) 
5  Maximum  ka  (integer) 
8  NREGNS  (integer) 


CYLFLPT . PRM 
1024 
64 


INPUTDAT . PRM 

•  ID-1  DKAO ,  minimum  ka  (REAL)  in  the  computation 

.875D-2  DINCA,  increment  of  ka  (REAL) 

40  NSTR,  start  point  (INTEGER)  in  the  computation 

479  NEND,  end  point  (INTEGER)  in  the  computation 

6. DO  RHA,  ratio  of  h  and  a-  (REAL) 

1  IE,  if  incident  wave  is  TE-polarized  it  is  1,  otherwise  0 

0  IM,  if  incident  wave  is  TM-polarized  it  is  1,  otherwise  0 

0-D0  THETAI ,  incident  angle  (REAL)  is  limited  to  0  to  90  degree 

6  NTHTAO ,  total  number  (INTEGER)  of  angle  of  THATA 

3  NPHI ,  total  number  (INTEGER)  of  angle  of  PHI0  to  be  computed 

0.D0  THETAO ,  initial  theta  angle  (REAL) 

1-D0  DELTHO,  increment  of  theta  angle  (REAL) 

0.DO  PHIA,  initial  theta  angle  (REAL) 

12. DO  DELPHI,  increment  of  phi  angle  (REAL) 

0  IK,  set  1  to  compute  scattering  currents  on  outer  and  inner  surface 


floating  point  zero  I0BIT  (1024  bits) 
floating  point  precision  IFPBIT  (64  bits) 


IMPEDNCE . PRM 

0  IZ,  if  perfect  conducting,  IZ=0;  otherwise,  IZ=1. 

( . 5D0 , 0 . DO )  Impedance  Z+  (phi  phi) 

( . 4D0 , 0 . DO )  Impedance  Z-  (phi  phi) 

( . 3D0 , 0 . DO )  Impedance  Z+  (phi  z) 

( . 3D0 , 0 . DO )  Impedance  Z-  (phi  z) 

(.3D0,0.D0)  Impedance  Z+  (z  phi) 
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(  . 3D0 , 0 . DO ) 
(  .4D0,0.D0) 
(  . 5D0 , 0 . DO ) 


Impedance  Z-  (z  phi) 
Impedance  Z+  (z  z) 
Impedance  Z-  (z  z) 


C.  SETUP  PROGRAM  AND  CREATED  FILES 

PROGRAM  SETUP 

C* ********************************************************************* 

C  NOTES: The  format  statement  1001  need  to  be  revised  if  other  than 
C  double  precision  real  numbers  are  used  for  DKHMAX  and  DKAMAX . 

C  Statement  1001  needs  to  be  revised  if  KHMAX  or  KAMAX  exceeds 

C  3  digits. 

C  The  format  statement  1002  need  to  be  revised  if  other  than 

C  double  precision  real  numbers  are  used  for  FZERO  AND  PRECSN. 

C  Statement  1002  needs  to  be  revised  if  I0BIT  exceeds  6  digits 

C  or  IFPBIT  exceeds  4  digits. 

C* ********************************************************************* 

c 

INCLUDE  ' REALTP . INC ' 

INCLUDE  ' CMPXTP . INC  1 

OPEN  (20, FILE= ' CYLGEOM . PRM 1 , IOSTAT=IOS , STATUS = ' OLD  1 ) 

IF ( IOS  .NE.  0)  THEN 
WRITE(*, 9000) 

9000  FORMAT  {'Cannot  find  the  file  CYLGEOM. PRM  containing  the  ',/, 

+  'maximum  values  for  ka  and  kh,  and  the  parameter  NREGNS . ') 

STOP 

END  IF 

READ  (20,*)  KHMAX 
READ  (20,*)  KAMAX 
READ  (20,*)  NREGNS 
CLOSE  (20) 

OPEN  (20, FILE= ' CYLFLPT . PRM ' , I OSTAT= IOS , STATUS =' OLD ' ) 

IF (IOS  .NE.  0)  THEN 
WRITE(*, 9001) 

9001  FORMAT  ('Cannot  find  the  file  CYLFLPT. PRM  containing  floating' ,/, 

+  ’  point  zero  bit  I0BIT  and  precision  IFPBIT.’) 

STOP 

END  IF 

READ (20,*)  I0BIT 
READ (20,*)  IFPBIT 
CLOSE  (20) 

OPEN  ( 2 1,FILE=' LIMITS. INC’ , STATUS =' UNKNOWN ' ) 

WRITE  (21,1001)  KHMAX,  KAMAX 

WRITE  (21,1002)  I0BIT,  IFPBIT,  IFPBIT 

CLOSE  (21) 

1001  FORMAT  (6X, 'PARAMETER  ( DKHMAX =  ’,13,’. DO,  DKAMAX=  ',13,'. DO)') 

1002  FORMAT  ( 6X ,' PARAMETER  (FZERO= ’ , 16 , ' . DO ,  PRECSN= ' , 14 , ' . DO , ' , 

+  '  IFPBIT= 1 ,14, ' ) ') 

C 

C  Part  2 

DKHMAX=KHMAX 
DKAMAX=KAMAX 
ONEDEL =ONE- EPS 8 

KQDIM=INT ( (DKHMAX /PI ) * NREGN S + ONEDEL ) 

KNDIM=INT( (DKAMAX/TWO) *NREGNS+ ONEDEL) 

KQDIM=MAX ( KQDIM , 1 ) 

KNDIM=MAX  ( KNDIM ,  1 ) 

OPEN  (21, FILE= ' MAINDM . INC ' , STATUS =' UNKNOWN ' ) 

WRITE  (21,1003)  NREGNS,  KNDIM,  KQDIM 
WRITE  (21,1004) 

WRITE  (21,1005) 

WRIE  (21,1006) 

CLOSE  (21) 

1003  FORMAT  ( 6X, ' PARAMETER  (NREGNS= ’ , 12 , ' ,  KNDIM= 1 , 13 , 

+  ’ ,  KQDIM= ' , 13 , 1 ) ’ ) 
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1004  FORMAT  ( 6X,  1  PARAMETER  ( KNDIM1 =KNDIM+ 1 ,  KQDIM1=KQDIM+1 )  ' ) 

1005  FORMAT  ( 6X, ' PARAMETER  (KXCRT=2 *KQDIMl ,  KCRNT=4 *KQDIMl ) 1 ) 

1006  FORMAT  ( 6X , ' PARAMETER  (MAXNG=KNDIM1 ,  MAXPEG=KQDIM/2+l , ' , 

+  1  MAXPOG=KQDIMl/2 ) ' ) 

C 

C  Part  3 

FZERO=I0BIT 

REFC=F ZERO* LOG (TWO) “LOG (PI2 ) -ONE 
REFH=HALF* {REFC-LOG (DKHMAX) ) 

REFA=HALF* (REFC-LOG (DKAMAX) ) 

C 

KQ2=KQDIM+2 

DMXM=KQ2 

DO  100  WHILE  ( (DMXM+HALF) * (LOG (DMXM/ DKHMAX) +ONEN)  . LT .  REFH) 
DMXM=  DMXM+ ONE 
100  CONTINUE 

MXMREG = INT { DMXM ) 

C 

MXMSNG =INT { TWO* DKHMAX +ONEDEL ) 

MXMSNG=MAX ( IFPBIT , MXMSNG , KQ2 ) +KNDIM+1 
OPEN  ( 21 , FILE= 1 GPQNDM . INC  1 , STATU S- ' UNKNOWN  1 ) 

WRITE  (21,1009)  MXMREG,  MXMSNG 
CLOSE  (21) 

1009  FORMAT  ( 6 X, ‘ PARAMETER  (MXMREG= ',14 , ' ,  MXMSNG= ' , 14 , ' ) ’ ) 

STOP 

END 


GPQNDM . INC 

PARAMETER  ( DKHMAX =  30. DO,  DKAMAX =  5. DO) 

PARAMETER  (FZERO=  1024. DO,  PRECSN=  64. DO,  IFPBIT=  64) 


LIMITS . INC 

PARAMETER  ( DKHMAX =  30. DO,  DKAMAX =  5. DO) 

PARAMETER  (FZERO=  1024. DO,  PRECSN=  64. DO,  IFPBIT=  64) 


MAINDM . INC 

PARAMETER  (NREGNS=  8,  KNDIM=  20,  KQDIM=  77) 

PARAMETER  (KNDIM1=KNDIM+1 ,  KQDIM1=KQDIM+1 ) 

PARAMETER  (KXCRT=2*KQDIM1 ,  KCRNT=4*KQDIM1 ) 

PARAMETER  (MAXNG=KNDIM1 ,  MAXPEG=KQDIM/2+l ,  MAXPOG=KQDIMl/2 ) 


D.  PROGRAM  MAIN 

PROGRAM  MAIN 

********************************************************************* 
INCLUDE  ' REALTP . INC ' 

INCLUDE  ' CMPXTP . INC ' 

COMMON  / GCONST/  DKH, DKA, HH, HA, HSQ, DASQ, HSQN, DASQN, HHSQ, HDASQ, 

+  RAH, RAHSQ, DHA, DAH 

COMMON  / INPUT 1/  DKA0 , DINCA, NSTR, NEND, RHA 
COMMON  / INPUT2 /  IE, IM, THETAI , THESINI , THECOSI , RTHEI 
COMMON  / INPUT4 /  IZ , IK, IS , NYSM 
C 

CALL  CHKINPUT 

OPEN  (21, FILE= ' rhz  z4 1 . dz  * , STATUS^ ' UNKNOWN ' ) 

DO  8001  I S  =NSTR , NEND 
DS=IS 

DKA=DKA0+ DINCA* DS 
DKH = DKA  *  RHA 
CALL  MAXODM 
CALL  RAPQMT 
CALL  GNPQFN 

CALL  BCKSCFL { IR , CESTH , CESPH ) 
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CALL  RCSPAREA (CESTH, CESPH) 
8001  CONTINUE 
CLOSE (21) 

STOP 

END 


E.  SUBROUTINE  CHKINPUT 


SUBROUTINE  CHKINPUT 

C******************************************************************************** 

INCLUDE  ' RE ALT P . INC ' 

INCLUDE  'LIMITS. INC' 

COMMON  / INPUT1 /  DKAO , DINCA  ,  NSTR , NEND , RHA 

COMMON  /INPUT2 /  IE, IM, THETAI , THESINI , THECOSI , RTHEI 

COMMON  /INPUT 3 /  RTHE , RDELT , RPHI , RDELP , NTHTAO , NPHI , THESIN , THECOS , 

+  RHP  I 

COMMON  /INPUT 4/  IZ, IK, IS,NYSM 

£*********************************************************************** 

OPEN (20, FILE= ' INPUTDAT . PRM '  , IOSTAT=IOS , STATUS= ' OLD  1 ) 

IF(IOS  -NE.  0)  THEN 

WRITE (*,*)  'Fail  to  open  input  file  INPUTDAT . PRM ’ 

STOP 

END  IF 

Q*»********************************************************************* 

C  Input  kh  and  ka.  These  values  are  passed  to  other 
C  parts  of  this  program  through  the  common  block  /INPUT1/. 

READ (20,*)  DKAO 
READ (20,*)  DINCA 
READ (20,*)  NSTR 
READ (20,*)  NEND 
READ (20,*)  RHA 

C****** ****************************************************** *********** 

C  Check  against  maximum  kh  and  ka  values. 

NBTW=NEND-NSTR 

DKH0=DKA0*RHA 

DKA1=DINCA*NBTW+DKA0 

DKH1=DKA1*RHA 

IF  (DKH1  .GT.  DKHMAX)  THEN 

IF  (DKA1  .GT.  DKAMAX)  THEN 

WRITE (*,*)  'Both  kh  and  ka  values  exceed  the  maximum  allowed.' 

ELSE 

WRITE (*,*)  'The  input  kh  value  exceeds  the  maximum  value  allowed.' 

END  IF 

WRITE (*,*)  'The  execution  is  terminated.' 

CLOSE (20) 

STOP 

ELSE  IF  (DKA1  . GT .  DKAMAX) THEN 

WRITE (*,*)  'The  input  ka  value  exceeds  the  maximum  value  allowed. ' 

WRITE ( * , * )  'The  execution  is  terminated.' 

CLOSE (20) 

STOP 

END  IF 

IF  ( (DINCH  .LT.  ZERO)  .OR.  (DINCA  .LT.  ZERO))  THEN 
WRITE  (*,*)  'The  increment  DINCH  or  DINCA  is  less  than  zero.' 

WRITE ( * , * )  'The  execution  is  terminated.' 

CLOSE (20) 

STOP 

END  IF 

£*********************************************************************** 

C  Input  the  incident  angle  (THETAI)  and  polarization  (TE  or  TM)  of  the 
C  incident  wave.  The  incident  angle  is  limited  to  0  to  90 
C  degrees.  SIN (THETAI)  and  COS (THETAI)  are  computed  also.  These 
C  parameters  are  passed  to  other  parts  of  this  program  through  the 
C  common  block  /INPUT2 / . 

READ (20,*)  IE 
READ (20,*)  IM 
READ (20,*)  THETAI 
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C  Check  the  input  value  of  incident  wave 

IF  ((IE  .NE.  1)  .AND.  (IE.  NE.  0))  THEN 
WRITE  (*,*)  'Improperly  specified  the  polarization  of  incident  ', 
+  'wave,  program  is  stoped. ' 

STOP 

ELSE  IF  ((IM  .NE.  1)  .AND.  (IM  .NE.  0))  THEN 
WRITE  {*,*)  'Improperly  specified  the  polarization  of  incident 
+  'wave,  program  is  stopped. ' 

CLOSE (20) 

STOP 

END  IF 


IF  ( (THETAI  .LT.  ZERO)  .OR.  (THETAI  .GT.  90.))  THEN 
WRITE  (*,*)  'Improperly  specified  incident  angle.  ', 

+  'program  is  stopped. ' 

CLOSE (20) 

STOP 

END  IF 

C  Calculate  SIN (THETAI)  and  COS (THETAI) 

RTHEI =THETAI *  PI / DEGPI 
THESINI=SIN ( RTHEI ) 

THECOSI =COS ( RTHEI ) 

£*★*** ****************************************** ************************ 

C  Input  the  angles  theta  and  phi  at  which  the  scattered  fields  are 
C  to  be  computed.  They  are  specified  in  terms  of  the  initial  theta 
C  (THETAO)  and  phi  (PHIO)  angles,  their  respective  increments  DELTHO 
C  and  DELPHI,  and  the  total  numbers  of  angles  NTHTAO  and  NPHI  to  be 
C  computed.  Thus  NTHTAO  and  NPHI  must  be  integers  greater  than  1.  If 
C  either  NTHTAO=0  or  NPHI=0,  no  bistatically  scattered  fields  will  be 
C  computed.  Note  that  the  scattered  electric  field  components  are 
C  computed  for  all  the  phi-angles  at  a  fixed  theta,  before  the 
C  theta-angle  is  varied.  All  angles  are  specified  in  degrees.  Theta  is 
C  limited  to  0  to  180  while  phi  is  limited  to  0  to  360  degrees. 

C  SIN (THETAO)  and  COS (THETAO)  are  computed  also.  These  parameters  are 
C  passed  to  other  parts  of  this  program  through  the  common  block  /INPUT3 / . 
READ (20,*)  NTHTAO 
READ (20,*)  NPHI 

IF  ((NTHTAO  .LT.  0)  .OR.  (NPHI  .LT.  0))  THEN 
WRITE (*,*)  'Improperly  specified  number  of  output  angles. 

+  'Program  is  stopped.' 

STOP 

ELSE  IF  ((NTHTAO  .EQ.  0)  .OR.  (NPHI  .EQ.  0))  THEN 
CLOSE (20) 

WRITE (*,*)  'Desired  bistatic  scattered  field  direction  has  not  ', 

+  'been  (properly)  specified,  they  will  not  be  computed. ' 

NTHTAO =0 
NPHI=0 
THETAO = ZERO 
DELTHO=ZERO 
PHIO=ZERO 
DELPHI = ZERO 
ELSE 

READ (20,*)  THETAO 
READ (20,*)  DELTHO 
READ (20, * )  PHIA 
READ (20,*)  DELPHI 
END  IF 

READ (20, * )  IK 


C  Check  input  values. 

IF  { (DKH0  .LE.  ZERO) 
WRITE (*,*)  'Invalid  kh, 
STOP 

END  IF 

IF  ( (DKA0  .LE.  ZERO) 
WRITE (*,*)  'Invalid  ka, 
STOP 

END  IF 


.OR.  (DKH1  .LE.  ZERO)) THEN 
program  is  stopped. ' 


.OR.  (DKA1  .LE.  ZERO ) ) THEN 
program  is  stopped. ' 
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IF  ( (DKAO  .GT.  DKHO)  .OR.  (DKA1  .GT.  DKH1 ) ) THEN 
WRITE{*,*)  'ka/kh  >  1,  program  is  stopped.' 

STOP 

END  IF 

C  Output  angle  checking  not  required: 

IF  ( NTHTAO  .EQ.  0)  GO  TO  200 
C  Chcking  output  angles : 

IF  ( (THETAO  .LT.  ZERO)  .OR.  (THETAO  .GT.  DEGPI) )  THEN 
WRITE (*,*)  'The  first  output  theta-angle  lies  outside  the  0  to 
+  '180  degrees  range.  Program  is  stopped.' 

STOP 

END  IF 

THTAIF=THETAO+ ( NTHTAO- 1 ) *DELTHO 

IF  { (THTAIF  .LT.  ZERO)  .OR.  (THTAIF  . GT .  DEGPI))  THEN 
WRITE (*,*)  'Some  of  the  specified  output  theta-angles  lie', 

+  'outside  the  0  to  180  degrees  range.  Program  is  stopped.' 
STOP 

END  IF 

PHIMX=TWO*DEGPI 

IF  ((PHIO  .LT.  ZERO)  .OR.  (PHIO  .GT.  PHIMX) )  THEN 
WRITE (*,*)  'The  first  output  phi-angle  lies  outside  the  0  to  ' 
+  '360  degrees  range.  Program  is  stopped.' 

STOP 

END  IF 

PHIF=PHIO+ (NPHI-1) * DEL PH I 

IF  ( (PHIF  .LT.  ZERO)  .OR.  (PHIF  . GT .  TWO*DEGPI))  THEN 
WRITE(*,*)  'Some  of  the  output  phi-angles  lie  outside', 

+  'the  0  to  360  degrees  range.  Program  is  stopped.' 

STOP 

END  IF 

200  CONTINUE 

DPI =PI /DEGPI 

RTHE=THETAO*DPI 

RPHI=PHIO*DPI 

RDELT=DELTHO*DPI 

RDELP=DELPHI*DPI 

RHPI=90.*DPI 

THESIN=SIN (RTHE) 

THECOS=COS (RTHE) 

RETURN 

END 


F.  SUBROUTINE  MAXODM 


SUBROUTINE  MAXODM 

Q  ********************************************************************* 

INCLUDE  ’ REALTP . INC ' 

INCLUDE  ' MAINDM . INC ' 

COMMON  /CRNTDM/  NMAX ,  MXNG ,  IQMAX ,  IQMAX1 ,  IQMAX2  ,  IXCRNT ,  ICRNT ,  MXQEG , 
+  MXQOG 

COMMON  /GCONST/  DKH , DKA , HH , HA , HSQ , DASQ , HSQN, DASQN , HHSQ , HDASQ , 

+  RAH , RAHSQ , DHA , DAH 

COMMON  /RTHETA/  DL1COSI , DL2SINI , DLL 
COMMON  /INPUT1 /  DKAO , DINCA, NSTR , NEND, RHA 
COMMON  /INPUT 2/  IE, IM, THETAI , THESINI , THECOSI , RTHEI 
COMMON  /INPUT3 /  RTHE , RDELT , RPHI , RDELP, NTHTAO , NPHI , THES IN, THECOS , 

+  RHPI 

C 

ONEDEL=ONE-EPS8 

IQMAX=INT ( (DKH/PI) *NREGNS+ONEDEL) 

IQMAX=MAX ( IQMAX , 1 ) 

IQMAX1=IQMAX+1 
IQMAX2=IQMAXl+2 
IXCRNT=2*IQMAXl 
ICRNT=4*IQMAX1 
MXQEG= I QMAX / 2  + 1 
MXQOG=IQMAXl / 2 
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c 

NMAX=INT ( (DKA/TWO) *NREGNS+ONEDEL) 

NMAX=MAX (NMAX, 1 ) 

MXNG=NMAX+ 1 
C 

C  Evaluate  SIN  and  COS  functions  to  pass  along  / 

DL1C0SI=DKH*THEC0SI 

DL2SINI=DKA*THESINI 

DLL=DKH*DL2SINI 

C 

C  Evaluate  geometrical  values  to  pass  along  /GCONST/ . 

HH=HALF*DKH 
HA=HALF  *  DKA 
HSQ=DKH*DKH 
D  AS  Q=  DKA  *  DKA 
HSQN=-HSQ 
DASQN=-DASQ 
HHSQ=QUAR  *HSQ 
HDASQ=QUAR*DASQ 
RAH=DKA/DKH 
RAHSQ=RAH*RAH 
DHA=DKH*HA 
DAH=HH*HA 
C 

RETURN 

END 

G.  SUBROUTINE  RAPQMT 

SUBROUTINE  RAPQMT 

C* ********************************************************************* 

INCLUDE  ' REALTP . INC ' 

INCLUDE  ' CMPXTP . INC ' 

INCLUDE  ' MAINDM . INC ' 

DIMENSION  00(2/2),  CIR(2,2) 

DIMENSION  CZSUM (2,2) ,CZDIF(2,2) ,  CR1 (4,4) , CR2 (4,4) , CR3 (4,4) 

DIMENSION  CRPQ1 ( KCRNT , KCRNT ) , CRPQ2 ( KCRNT , KCRNT ) , 

+  CRPQ3 1 ( KCRNT , KCRNT ) , CRPQ3 2 ( KCRNT , KCRNT ) 

COMMON  /INPUT 2/  IE, IM, THETAI , THESINI , THECOSI , RTHEI 
COMMON  /INPUT 4/  IZ , IK, IS , NYSM 
COMMON  /XPQTMP1 /  CRPQ1 , CRPQ2 
COMMON  /XPQTMP2 /  CRPQ31 , CRPQ32 

COMMON  /CRNTDM/  NMAX, MXNG , IQMAX, IQMAX1 , IQMAX2 , IXCRNT, ICRNT , MXQEG, 

+  MXQOG 

c 

OPEN  (20 , FILE= ' IMPEDNCE . PRM ' , IOSTAT=IOS , STATUS- 1  OLD 1 ) 

IF  (IOS  .NE.  0)  THEN 
WRITE ( * , 9000) 

9000  FORMAT  ('Cannot  find  the  file  IMPEDNCE 1 . PRM  containing  the  ’,/, 

+  ' impedances  of  the  inner  and  outer  surfaces . ' ) 

STOP 

END  IF 
C 

READ (20,*)  IZ 

IF  {(IZ  .NE.  0)  .AND.  (IZ  .NE.  1))  THEN 
WRITE  (*,*)  'Improperly  specified  the  impedance  of  cylinder  ', 

+  'surface,  program  is  stopped.1 
CLOSE (20) 

STOP 

END  IF 

IF  (IZ  .EQ.  0)  THEN 
CLOSE (20) 

RETURN 
END  IF 

C  Set  up  the  matrix  when  Z  and  Delta  are  diagonal,  or  not  diagonal  but  n  <0. 
C  Initialize  the  matrix  CRPQ1 
DO  200  1=1, KCRNT 
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DO  100  J=1 , KCRNT 
CRPQ1 ( I , J ) =CZERO 
CRPQ2(I, J)=CZERO 
CRPQ31 (I, J)=CZERO 
CRPQ32 (I, J)=CZERO 
100  CONTINUE 

200  CONTINUE 

DO  400  1=1,2 
DO  300  J=l,2 
READ  (20,*)  CO(I,J) 

READ  (20,*)  CIR ( I , J) 

CZSUM { I , J) = (CO ( I ,  J) +CIR ( I , J) ) *HALF 
CZDIF < I , J)  =  (CO ( I ,  J) -CIR ( I , J) ) *HALF 
300  CONTINUE 

400  CONTINUE 

CLOSE  (20) 

CZ11=CZSUM(1,1) 

CZ12=CZSUM (1,2) 

CZ21=CZSUM (2,1) 

CZ22=CZSUM (2,2) 

CRDET=CZ11*CZ22-CZ12*CZ21 
CD11=CZDIF (1,1) 

CD12=CZDIF (1,2) 

CD2 1=CZDIF (2,1) 

CD22=CZDIF (2,2) 

IF  ( CRDET  .EQ.  CZERO)  THEN 
WRITE  (*,9001) 

STOP 

END  IF 

C  Check  the  diagonalization  of  the  impedance  matrixes 
NSYM=1 

IF  (CZ12  .NE .  CZERO)  THEN 
NSYM=0 

ELSE  IF  (CZ21  .NE.  CZERO)  THEN 
NSYM=0 

ELSE  IF  (CD12  .NE.  CZERO)  THEN 
NSYM=0 

ELSE  IF  (CD21  .NE.  CZERO)  THEN 
NSYM=0 
END  IF 

CRDT Z = CONE / CRDET 

CR1 (1, 1) =CZ11- (CD11*CD11*CZ22-CD11*CD12*CZ21-CD11*CD21*CZ12 
+  +CD12 *CD21*CZ11 ) *CRDTZ 

CR1 (1, 2)=CZ12- (CD11*CD12*CZ22-CD12*CD12*CZ21-CD11*CD22*CZ12 
+  +CD12*CD22*CZ11) *CRDTZ 

CR1 (2,1) =CZ21- (CD11*CD21*CZ22-CD11*CD22*CZ21-CD21*CD21*CZ12 
+  +CD21*CD22*CZ11)*CRDTZ 

CR1 (2,2) =CZ22- (CD12  *CD21*CZ22-CD12*CD22*CZ21-CD21*CD22*CZ12 
+  +CD22*CD22*CZ11) *CRDTZ 

CR1 (1,3)= (CD12*CZ11-CD11*CZ12 ) *CRDTZ 
CR1 (1, 4)= (CD12*CZ21-CD11*CZ22) *CRDTZ 
CR1 (2,3)= (CD22 *CZ11-CD21*CZ12 ) *CRDTZ 
CR1 (2, 4) = (CD22*CZ21-CD21*CZ22) *CRDTZ 
CR1 (3,1)= (CD11*CZ21-CD21*CZ11 ) *CRDTZ 
CR1 (3,2)= (CD12*CZ21-CD22*CZ11 ) *CRDTZ 
CR1 (4,1)= (CD11*CZ22 -CD21*CZ12 ) *CRDTZ 
CRl (4,2)= (CD12*CZ22-CD22*CZ12 ) *CRDTZ 
CR1 (3,3) =CZ11*CRDTZ 
CRl (3,4) =CZ21*CRDTZ 
CRl (4,3) =CZ12*CRDTZ 
CRl (4,4) =CZ22*CRDTZ 
C 

DO  1300  IQE=0 , MXQEG 

IQE1=IQE+1 

IQ=2*IQE 

DQ=IQ 

DQ1=IQ+1 

IQX1=4*IQ+1 
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IQX2=IQX1+1 

IQX3=IQX2+1 

IQX4=IQX3+1 

DO  1100  IPE=0 , MXQEG 
IPEl=IPE+l 
IP=2*IPE 
IPX1=4*IP+1 
IPX2=IPX1+1 
IPX3=IPX2+1 
IPX4=IPX3+1 
DP=IP 
DP1=IP+1 

DBPHI=PISQ* (DP+DQ1 ) * (DP1-DQ) 

BPHI=4 .DO/DBPHI 

DBZ=PISQ* (DP+DQ1) * (DP1-DQ) * (DP1+DQ1+ONE ) * (DP-DQ1) 
BZ  =  -*8  .  D0*DQ1/DBZ 
CRPQ1 ( IPX1 , IQX1 ) =CR1 (1,1) *BPHI 
CRPQ1 ( IPX2 , IQX1 ) =CR1 (2,1) *BPHI 
CRPQ1 ( IPX3 , IQX1 ) =CR1 (3,1) *BPHI 
CRPQ1 ( IPX4 , IQX1 ) =CR1 (4,1) *BPHI 
CRPQ1 ( IPX1 , IQX2 ) =CR1 ( 1 , 2 ) *BZ 
CRPQ1 (IPX2 , IQX2 ) =CR1 (2 , 2 ) *BZ 
CRPQ1 ( IPX3 , IQX2 ) =CR1 ( 3 , 2 ) *BZ 
CRPQ1 ( IPX 4 , IQX2 ) =CR1 ( 4 , 2 ) *BZ 
CRPQl ( IPX1 , IQX3 ) =CR1 (1,3) *BPHI 
CRPQ1 (IPX2 , IQX3 ) =CR1 (2,3) *BPHI 
CRPQl ( IPX3 , IQX3 ) =CR1 (3,3) *BPHI 
CRPQl ( IPX4 , IQX3 ) =CR1 (4,3) *BPHI 
CRPQl ( IPX1 , IQX4 ) =CR1 (1,4) *BZ 
CRPQl ( IPX2 , IQX4 ) =CR1 (2,4) *BZ 
CRPQl ( IPX3 , IQX4 ) =CR1 ( 3 , 4 ) *BZ 
CRPQl ( IPX4 , IQX4 ) =CR1 ( 4 , 4 ) *BZ 
1100  CONTINUE 

1300  CONTINUE 

DO  2300  IQO=0 , MXQOG 

IQ01=IQ0+1 

IQ=2*IQ0+1 

IQX1=4*IQ+1 

IQX2=IQX1+1 

IQX3=IQX2+1 

IQX4=IQX3+1 

DQ=IQ 

DQ1=IQ+1 

DO  2200  IPO=0, MXQOG 
IP01=IPO+l 
IP=2*IPO+l 
IPX1=4*IP+1 
IPX2=IPX1+1 
IPX3=IPX2+1 
IPX4=IPX3+1 
DP-IP 
DP1=IP+1 

DBPHI=PISQ* (DP+DQ1) * (DP1-DQ) 

BPHI= 4. DO/DBPHI 

DBZ=PISQ* (DP+DQ1 ) * (DPl-DQ) * (DP1+DQ1+0NE) * (DP-DQ1) 

BZ=-8 . D0*DQ1 /DBZ 

CRPQl ( IPX1 , IQX1 ) =CR1 (1,1) *BPHI 

CRPQl ( IPX2 , IQX1 ) =CRl (2,1) *BPHI 

CRPQl ( IPX3 , IQX1 ) =CR1 (3,1) *BPHI 

CRPQl ( IPX4 , IQX1 ) =CR1 (4,1) *BPHI 

CRPQl ( IPX1 , IQX2 ) =CR1 ( 1 , 2 ) *BZ 

CRPQl ( IPX2 , IQX2 ) =CR1 ( 2 , 2 ) *BZ 

CRPQl ( IPX3 , IQX2 ) =CR1 { 3 , 2 ) *BZ 

CRPQl ( IPX4 , IQX2 ) =CRl ( 4 , 2 ) *BZ 

CRPQl ( IPX1 , IQX3 ) =CR1 (1,3) *BPHI 

CRPQl ( IPX2 , IQX3 ) =CR1 (2,3) *BPHI 

CRPQl ( IPX3 , IQX3 ) =CR1 (3,3) *BPHI 

CRPQl ( IPX4 , IQX3 ) =CR1 (4,3) *BPHI 
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CRPQ1 (IPX1 , IQX4 ) =CR1 (1 , 4) *BZ 
CRPQ1 ( IPX2 , IQX4 ) =CRl ( 2 , 4 ) *BZ 
CRPQ1 ( IPX3 , IQX4 ) =CR1 ( 3 , 4 ) *  BZ 
CRPQ1 (IPX4 , IQX4 ) =CR1 {4 , 4 ) *BZ 
2200  CONTINUE 

2300  CONTINUE 

C  Set  up  the  matrix  Rpq  utilized  when  Z  or  Delta  is  nor  diagonal,  and  n  <  0 
CR2 (1,2) =-CRl (1,2) 

CR2 (2,1) =-CRl (2,1) 

CR2 (1,3) =  -CRl (1,3) 

CR2 (2,4) =-CRl (2,4) 

CR2 (3,1) =-CRl (3,1) 

CR2 (4,2) =-CRl (4,2) 

CR2 (3,4)= -CR1 (3,4) 

CR2 (4,3) =-CRl (4,3) 

CR2 (1,1) =CR1 (1,1) 

CR2 (2,2) =CR1 (2,2) 

CR2 (1,4) =CR1 (1,4) 

CR2 (2,3) =CR1 (2,3) 

CR2 (3,2) =CR1 (3,2) 

CR2 (4,1) =CR1 (4,1) 

CR2 (3,3) =CR1 (3,3) 

CR2 (4,4) =CR1 (4,4) 

DO  3300  IQE=0 , MXQEG 

IQE1=IQE+1 

IQ=2*IQE 

DQ=IQ 

DQ1=IQ+1 

IQX1=4*IQ+1 

IQX2=IQX1+1 

IQX3=IQX2+1 

IQX4=IQX3+1 

DO  3100  IPE=0, MXQEG 
IPE1=IPE+1 
IP=2  *IPE 
IPXl  =  4  *IP+1 
IPX2=IPX1+1 
IPX3=IPX2+1 
IPX4=IPX3  +1 
DP=IP 
DP1=IP+1 

DBPHI=PISQ* (DP+DQl) * (DP1-DQ) 

BPHI=4 . D0/DBPHI 

DBZ^PISQ* (DP+DQ1) * (DP1-DQ) * (DP1+DQ1+ONE) * (DP-DQ1) 

BZ= - 8 . DO  *DQ1 /DBZ 
CRPQ2 ( IPX1 , IQX1 ) =CR2 (1,1) *BPHI 
CRPQ2 ( IPX2 , IQXl ) =CR2 (2,1) *BPHI 
CRPQ2 ( IPX3 , IQXl ) =CR2 (3,1)*BPHI 
CRPQ2 ( IPX4 , IQXl ) =CR2 (4,1) *BPHI 
CRPQ2 ( IPX1 , IQX2 ) =CR2 ( 1 , 2 ) *BZ 
CRPQ2 ( IPX2 , IQX2 ) =CR2 ( 2 , 2 ) *BZ 
CRPQ2 { IPX3 , IQX2 ) =CR2 ( 3 , 2 ) *BZ 
CRPQ2 ( IPX4 , IQX2 ) =CR2 (4 , 2 ) *BZ 
CRPQ2 ( IPX1 , IQX3 ) =CR2 (1,3) *BPHI 
CRPQ2 (IPX2 , IQX3 ) =CR2 (2, 3) *BPHI 
CRPQ2 ( IPX3 , IQX3 ) =CR2 (3,3) *BPHI 
CRPQ2 ( IPX4 , IQX3 ) =CR2 (4,3) *BPHI 
CRPQ2 (IPX1, IQX4)=CR2 (1,4) *BZ 
CRPQ2 ( IPX2 , IQX4 ) =CR2 ( 2 , 4 ) *BZ 
CRPQ2 { IPX3 , IQX4 ) =CR2 ( 3 , 4 ) *BZ 
CRPQ2 { IPX4 , IQX4 ) =CR2 (4,4)*BZ 
3100  CONTINUE 

3300  CONTINUE 

DO  4300  IQO=0 , MXQOG 

IQO!=IQO+l 

IQ=2*IQO+l 

IQX1=4*IQ+1 

IQX2=IQX1+1 
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IQX3=IQX2+1 

IQX4=IQX3+1 

DQ=IQ 

DQ1=IQ+1 

DO  4200  IPO=0 , MXQOG 
IP01=IP0+1 
IP=2*IPO+l 
IPX1=4*IP+1 
IPX2=IPX1+1 
IPX3=IPX2+1 
IPX4=IPX3+1 
DP=IP 
DP1=IP+1 

DBPHI=PISQ* (DP+DQl ) * (DP1-DQ) 

BPHI=4 /DBPHI 

DBZ=PISQ* (DP+DQl ) * (DPl-DQ) * (DP1+DQ1+ONE) * (DP-DQl) 

BZ=-8.*DQ1/DBZ 

CRPQ2 ( IPX1 , IQX1 ) =CR2 (1,1)*BPHI 
CRPQ2 ( IPX2 , IQX1 ) =CR2 (2,1) *BPHI 
CRPQ2 ( IPX3 , IQX1 ) =CR2 (3,1) *BPHI 
CRPQ2 ( IPX4 , IQX1 ) =CR2 (4,1) *BPHI 
CRPQ2 ( IPX1 , IQX2 ) =CR2 ( 1 , 2 ) *BZ 
CRPQ2 ( IPX2 , IQX2 ) =CR2 ( 2 , 2 ) *BZ 
CRPQ2 (IPX3 , IQX2 ) =CR2 (3 , 2 ) *BZ 
CRPQ2 ( IPX 4 , IQX2 ) =CR2 ( 4 , 2 ) *BZ 
CRPQ2 ( IPX1 , IQX3 ) =CR2 (1,3) *BPHI 
CRPQ2 (IPX2 , IQX3 ) =CR2 (2,3) *BPHI 
CRPQ2 (IPX3 , IQX3 ) =CR2 (3,3) *BPHI 
CRPQ2 { IPX 4 , IQX3 ) =CR2 (4,3) *BPHI 
CRPQ2 { IPX1 , IQX4 ) =CR2 (1,4)*BZ 
CRPQ2 (IPX2 , XQX4) =CR2 (2 , 4) *BZ 
CRPQ2 (IPX3 , IQX4 ) =CR2 (3 , 4) *BZ 
CRPQ2 ( IPX4 , IQX4 ) =CR2 { 4 , 4 ) *BZ 
4200  CONTINUE 

4300  CONTINUE 

C  Prepare  a  mtarix  for  computing  the  scattering  currents  on  the  outer 
C  and  the  inner  surfaces 

IF  (IK  .NE.  1)  THEN 
RETURN 

END  IF 

CR3 (1,1) =-CRl (4,1) 

CR3 (1,2) =-CRl (4,2) 

CR3 (1,3) =-CRl (4,3) 

CR3 (1,4) =~CR1 (4,4) 

CR3 (3,1) =CR1 (2,1) 

CR3 (3,2) =CR1 (2,2) 

CR3 (3,3) =CR1 (2,3) 

CR3 (3,4) =CR1 (2,4) 

CR3 (2,1) =CR1 (3,1) 

CR3 (2,2) =CR1 (3,2) 

CR3 (2,3) =CR1 (3,3) 

CR3 (2,4) =CR1 (3,4) 

CR3 (4,1) =-CRl (1,1) 

CR3 (4,2) =“CR1 (1,2) 

CR3 (4,3) =-CRl (1,3) 

CR3 (4,4)= -CRl (1,4) 

C 

DO  5200  IQE=0 , MXQEG 

IQE1=IQE+1 

IQ=2*IQE 

DQ=IQ 

DQ1=IQ+1 

IQX1=4*IQ+1 

IQX2=IQX1+1 

IQX3=IQX2+ 1 

IQX4=IQX3+1 

DO  5100  IPE=0, MXQEG 
IPE1=IPE+1 
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IP=2  * IPE 

IPX1=4*IP+1 

IPX2=IPX1+1 

IPX3=IPX2+1 

IPX4=IPX3+1 

DP=IP 

DPl=IP+l 

DBPHI=PISQ* (DP+DQ1) * (DP1-DQ) 

BPHI=4 . /DBPHI 

DBZ=PISQ* (DP+DQ1) * (DPl-DQ) * (DP1+DQ1+0NE) * (DP-DQ1) 
BZ=-8 . *DQ1/DBZ 

CRPQ3 1 ( IPX1 , IQXl ) = { CONE+CR3 (1,1)) *BPHI 
CRPQ3 1 ( IPX2 , IQXl ) =CR3 (2,1) *BPHI 
CRPQ3 1 ( IPX3 , IQXl ) =CR3 (3,1) *BPHI 
CRPQ3 1 ( IPX4 , IQXl ) =CR3 (4,1) *BPHI 
CRPQ3 1 ( IPX1 , IQX2 ) =CR3 ( 1 , 2 ) *BZ 
CRPQ3 1 ( IPX2 , IQX2 ) = ( CONE+CR3 ( 2 , 2 ) ) *BZ 
CRPQ3 1 ( IPX3 , IQX2 ) =CR3 ( 3 , 2 ) *BZ 
CRPQ3 1 ( IPX4 , IQX2 ) =CR3 (4,2) *BZ 
CRPQ3 1 ( IPX1 , IQX3 ) =CR3 (1,3) *BPHI 
CRPQ3 1 ( IPX2 , IQX3 ) =CR3 (2,3) *BPHI 
CRPQ31 (IPX3 , IQX3 ) = (CONE+CR3 (3,3)) *BPHI 
CRPQ3 1 ( IPX4 , IQX3 ) =CR3 (4,3) *BPHI 
CRPQ3 1 ( IPX1 , IQX4 ) =CR3 ( 1 , 4 ) *BZ 
CRPQ31 (IPX2 , IQX4) =CR3 (2,4) *BZ 
CRPQ3 1 ( IPX3 , IQX4 ) =CR3 ( 3 , 4 ) *BZ 
CRPQ3 1 ( I PX4 , IQX4 ) = ( CONE+CR3 ( 4 , 4 ) ) *BZ 
C 

CRPQ3  2 ( IPXl , IQXl ) = ( CONE-CR3 (1,1)) *BPHI 
CRPQ32 (IPX2 , IQXl) =-CR3 (2,1) *BPHI 
CRPQ32 ( IPX3 , IQXl ) =-CR3 (3,1) *BPHI 
CRPQ32 ( IPX4 , IQXl ) =-CR3 (4,1)*BPHI 
CRPQ3 2 ( IPXl , IQX2 ) =-CR3 ( 1 , 2 ) *BZ 
CRPQ32 ( IPX2 , IQX2 ) = ( CONE-CR3 (2,2) )*BZ 
CRPQ32 ( IPX3 , IQX2 ) =-CR3 (3,2) *BZ 
CRPQ32 ( IPX4 , IQX2 ) =-CR3 (4,2) *BZ 
CRPQ32 (IPXl, IQX3)=-CR3 (1,3) *BPHI 
CRPQ32 ( IPX2 , IQX3 ) =-CR3 (2,3) *BPHI 
CRPQ32 (IPX3 , IQX3) = (CONE-CR3 (3,3)) *BPHI 
CRPQ32 ( IPX4 , IQX3 ) =-CR3 (4,3) *BPHI 
CRPQ3 2 ( IPXl , IQX4 ) =-CR3 ( 1 , 4 ) *BZ 
CRPQ32 ( IPX2 , IQX4 ) =-CR3 (2,4)*BZ 
CRPQ32 ( IPX3 , IQX4 ) =-CR3 (3 , 4) *BZ 
CRPQ32 ( IPX4 , IQX4 ) = ( CONE-CR3 (4,4) ) *BZ 
5100  CONTINUE 

5200  CONTINUE 
C 

DO  6200  IQO=0 , MXQOG 

IQ01=IQO+l 

IQ=2*IQO+l 

DQ=IQ 

DQ1=IQ+1 

IQX1=4*IQ+1 

IQX2=IQX1+1 

IQX3=IQX2+1 

IQX4=IQX3+1 

DO  6100  IPO=0, MXQOG 
IPEl=IPO+l 
IP=2*IPO+l 
IPX1=4*IP+1 
IPX2=IPX1+1 
IPX3=IPX2+1 
IPX4=IPX3+1 
DP=IP 
DP1=IP+1 

DBPHI=PISQ* (DP+DQ1 ) * (DPl-DQ) 

BPHI=4. /DBPHI 

DBZ=PISQ* (DP+DQ1) * (DPl-DQ) * (DP1+DQ1+ONE) * (DP-DQ1) 
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BZ=-8 . *DQ1/DBZ 

CRPQ31 (IPXl, IQX1)= (C0NE+CR3 (1,1)) *BPHI 
CRPQ3 1 ( IPX2 , IQX1 ) =CR3 (2,1) *BPHI 
CRPQ31 (IPX3 , IQX1) =CR3 (3 , 1) *BPHI 
CRPQ3 1 ( I PX4 , IQX1 ) =CR3 (4,1)*  BPHI 
CRPQ3 1 ( IPX1 , IQX2 ) =CR3 ( 1 , 2 ) *BZ 
CRPQ3 1 ( IPX2 , IQX2 ) = ( C0NE+CR3 ( 2 , 2 ) ) *BZ 
CRPQ3 1 ( IPX3 , IQX2 ) =CR3 ( 3 , 2 ) *BZ 
CRPQ3 1 ( IPX4 , IQX2 ) =CR3 ( 4 , 2 ) *BZ 
CRPQ3 1 { IPX1 , IQX3 ) =CR3 (1,3) *BPHI 
CRPQ3 1 ( IPX2 , IQX3 ) =CR3 (2,3) *BPHI 
CRPQ31 (IPX3 , IQX3 ) = (C0NE+CR3 (3,3)) *BPHI 
CRPQ3 1 ( IPX 4 , IQX3 ) =CR3 (4,3) *BPHI 
CRPQ3 1 ( IPX1 , IQX4 ) =CR3 ( 1 , 4 ) *BZ 
CRPQ3 1 ( IPX2 , IQX4 ) =CR3 ( 2 , 4 ) *BZ 
CRPQ31 ( IPX3 , IQX4 ) =CR3 ( 3 , 4 ) *BZ 
CRPQ3 1 ( IPX4 , IQX4 ) = ( C0NE+CR3 (4 , 4 ) ) *BZ 

C 

CRPQ32 ( IPXl , IQX1 ) = (C0NE-CR3 (1,1)) *BPHI 
CRPQ32 ( IPX2 , IQX1 ) =-CR3 (2,1) *BPHI 
CRPQ32 ( IPX3 , IQX1 ) =-CR3 ( 3 , 1 ) *BPHI 
CRPQ32 ( IPX4 , IQX1 ) =-CR3 ( 4 , 1 ) *BPHI 
CRPQ32 ( IPXl , IQX2 ) =-CR3 (1,2)*BZ 
CRPQ32 ( IPX2 , IQX2 ) = ( C0NE-CR3 (2,2) )*BZ 
CRPQ32 (IPX3 , IQX2 ) =-CR3 (3,2)*BZ 
CRPQ3 2 ( IPX4 , IQX2 ) = -CR3 ( 4 , 2 ) *BZ 
CRPQ32 (IPXl , IQX3 ) =-CR3 (1,3)*BPHI 
CRPQ32 ( IPX2 , IQX3 ) =-CR3 (2,3)*BPHI 
CRPQ32 (IPX3 , IQX3 ) = (C0NE-CR3 (3,3)) *BPHI 
CRPQ32 ( IPX4 , IQX3 ) =-CR3 (4,3) *BPHI 
CRPQ3 2 (IPXl , IQX4 ) = -CR3 ( 1 , 4 ) *BZ 
CRPQ32 (IPX2 , IQX4 ) =-CR3 (2,4) *BZ 
CRPQ3 2 ( IPX3 , IQX4 ) =-CR3 ( 3 , 4 ) *BZ 
CRPQ32 (IPX4 , IQX4 ) = (C0NE-CR3 (4,4) )*BZ 

6100  CONTINUE 

6200  CONTINUE 
RETURN 

9001  FORMAT (’The  given  inside  and  outside  impedance  have  a  singular', 
+  '  sum.  Execution  is  terminated. ' ) 

END 


H.  SUBROUTINE  GNPQFN 

SUBROUTINE  GNPQFN 

£*★******************★********************************************************** 
INCLUDE  *  REALTP . INC ' 

INCLUDE  '  MAINDM .  INC  1 

COMMON  /GCONST/  DKH, DKA, HH, HA, HSQ, DASQ, HSQN, DASQN, HHSQ, HDASQ, 

+  RAH , RAHSQ , DHA , DAH 

CHARACTER  FILEEVEN1*12 ,  FILEODDl*12 
C  Set  up  the  file  names 
NC1=INT (DKA) 

DC=NC1 

NC=100000 . *DKH+NC1 
DC1=DKA-DC 
NC2=INT (1000. *DC1 ) 

C  Set  up  the  indices  of  file  names  of  Green's  function,  and  check  whether 
C  these  files  exist.  If  these  files  exist,  then  use  it  directly.  Otherwise, 

C  call  subroutine  XPQINI1. 

FILEEVEN1='E 
FILE0DD1= ' O 

IGE=8*2* (MAXPEG+1 ) * (MAXPEG+2) 

IG0=8*2* (MAXPOG+l) * (MAXPOG+2) 

WRITE  (FILEEVEN1 (2:8), ' (17.7) ’ )  NC 
WRITE  (FILEODD1 (2:8) ,' (17.7) 1 )  NC 
WRITE  (FILEEVENl (10:12) , ' (13.3) ')  NC2 
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WRITE  (FILE0DD1 (10:12), ' (13.3) ')  NC2 

OPEN  (28, ACCESS^' DIRECT' , FILE=FILEEVEN1 , RECL=IGE , IOSTAT=IOS , 
+  STATUS= ' OLD  1 ) 

IF  (IOS  .NE.  0)  THEN 
CLOSE  (28) 

CALL  XPQINI 1 ( DKH , DKA ) 

ELSE 

OPEN  (29, ACCESS^ ' DIRECT  1  , FILE=FILEODDl , RECL=IGO, IOSTAT=IOS, 

+  STATUS= ' OLD  1 ) 

IF  (IOS  .NE.  0)  THEN 
CLOSE  (29) 

CALL  XPQINI1 (DKH, DKA) 

ELSE 
CLOSE  (28) 

CLOSE  (29) 

CALL  XPQINI (DKH, DKA) 

END  IF 
END  IF 
RETURN 
END 


C 

SUBROUTINE  XPQINI (DKHIN, DKAIN) 

Q* ********************************************************************** 

INCLUDE  1 REALTP . INC ' 

INCLUDE  ' CMPXTP . INC ' 

INCLUDE  ' MAINDM . INC  1 

DIMENSION  GNE ( 4 , (MAXPEG+1 ) * (MAXPEG+2 ) /2 ) , 

+  GNO (4 , (MAXPOG+1 ) * (MAXPOG+2 ) /2 ) 

DIMENSION  CGNE ( 0 : MAXPEG , 0 :  MAXPEG , KNDIM1  +  1 )  , 

+  CGNO ( 0 : MAXPOG , 0 : MAXPOG , KNDIMl + 1 ) 

DIMENSION  CDGNE ( 0 : MAXPEG , 0 : MAXPEG , KNDIMl + 1 ) , 

+  CDGNO ( 0 : MAXPOG , 0 : MAXPOG, KNDIMl +1) 

COMMON  / CRNTDM/  NMAX, MXNG, IQMAX, IQMAX1 , IQMAX2 , IXCRNT, ICRNT,MXQEG, 

+  MXQOG 

COMMON  / GPQTMP /  CGNE , CDGNE , CGNO , CDGNO 
SAVE  /GPQTMP/ 

C 

CHARACTER  FILEEVEN1*12 ,  FILEODDl*12 
DKH = DKH IN 
DKA= DKAIN 


C  Initialize  the  matrix  CGNE,  CGNO,  CDGNE,  CDGNO 
DO  105  IC=  1 , KNDIM1+1 

DO  102  IB=  0, MAXPEG 
DO  101  IA=  0 , MXPEG 
CGNE ( I A , IB , IC) =CZERO 
CDGNE (IA, IB, IC) =CZERO 

101  CONTINUE 

102  CONTINUE 

DO  104  ID=0, MAXPOG 
DO  103  IE=0, MAXPOG 
CGNO ( IE , ID , IC ) =CZERO 
CDGNO ( IE , ID , IC ) =CZERO 

103  CONTINUE 

104  CONTINUE 

105  CONTINUE 
C 

C  Set  up  the  file  name 
NC1=INT (DKA) 

DC=NC1 

NC=100000 . *DKH+NC1 
DC 1= DKA -DC 
NC2=INT(1000.*DC1) 

FILEEVEN1= ' E 
FILEODDl=  *  O 

IGE=8*2  * (MAXPEG+1) * (MAXPEG+2) 

IGO=8*2* (MAXPOG+1) * (MAXPOG+2) 

WRITE  (FILEEVEN1 (2 : 8)  ,  1  (17 .7)  '  )  NC 
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WRITE  ( FILEODDl (2:8), 1 (17.7) ')  NC 
WRITE  (FILEEVEN1 (10:12), ' (13.3) 1 )  NC2 
WRITE  (FILEODDl (10:12),  '  (13.3)  ' )  NC2 

OPEN  (28, ACCESS= 1  DIRECT ' , FILE=FILEEVEN1 , RECL=IGE, IOSTAT=IOS, 

+  STATU S= ' OLD ' ) 

C 

OPEN  (29,ACCESS= 'DIRECT' , FILE=FILEODDl , RECL=IGO, IOSTAT=IOS , 

+  STATUS= ' OLD ' ) 

DO  900  NI=1 , MXNG+1 

C  The  following  values  N,  DN  and  DNH  are  passed  to  the  G-computation 

C  related  subroutines  through  the  common  block  /NCONST/ . 

READ  (28 , REC=NI )  GNE 
IRECE=0 

DO  500  IQE=0 ,  MXQEG 
DO  300  IPE=0 , IQE-1 
CGNE ( I PE , IQE , NI ) =CGNE ( IQE , I PE , NI ) 

CDGNE ( IPE , IQE , NI ) =CDGNE ( IQE , IPE , NI ) 

300  CONTINUE 

DO  400  I PE = IQE, MXQEG 
IRECE= IRECE+ 1 
GR=GNE ( 1 , IRECE) 

GI=GNE (2 , IRECE) 

GDR=GNE (3 , IRECE) 

GDI=GNE(4, IRECE) 

CGNE ( IPE , IQE , NI ) =DCMPLX ( GR , GI ) 

CDGNE ( IPE , IQE , NI ) =DCMPLX (GDR , GDI ) 

400  CONTINUE 

500  CONTINUE 

READ  (29 , REC=NI )  GNO 
IRECO=Q 

DO  800  IQO=0 , MXQOG 
DO  600  IPO=0 , IQO-1 
CGNO ( I PO , IQO , NI ) =CGNO ( IQO , IPO , NI ) 

CDGNO ( IPO , IQO , NI ) =CDGNO ( IQO , IPO , NI ) 

600  CONTINUE 

DO  700  IPO=IQO, MXQOG 
IRECO=IRECO+l 
GR=GNO ( 1 , IRECO ) 

GI=GNO ( 2 , IRECO) 

GDR=GNO (3 , IRECO) 

GDI=GNO ( 4 , IRECO) 

CGNO ( I PO , IQO , NI ) =DCMPLX ( GR , GI ) 

CDGNO ( IPO , IQO , NI ) =DCMPLX ( GDR , GDI ) 

700  CONTINUE 

800  CONTINUE 

900  CONTINUE 
CLOSE  (28) 

CLOSE  (29) 

RETURN 

END 

C 

SUBROUTINE  XPQINIl (DKHIN, DRAIN) 

Q*********************************************************************** 

INCLUDE  *  RE ALT P . INC ' 

INCLUDE  ' CMPXTP . INC ' 

INCLUDE  ' MAINDM . INC ' 

INCLUDE  ’ GPQNDM . INC ' 

INCLUDE  'LIMITS. INC* 

DIMENSION  GNE (4, (MAXPEG+1 } * (MAXPEG+2 ) /2 ) , 

+  GNO (4, (MAXPOG+1 ) * (MAXPOG+2) /2) 

DIMENSION  CGNE (0 : MAXPEG, 0 : MAXPEG, KNDIMl +1)  , 

+  CGNO ( 0 : MAXPOG , 0 : MAXPOG , KNDIM1 +1 ) 

DIMENSION  CDGNE ( 0 : MAXPEG , 0 : MAXPEG , KNDIMl + 1 ) , 

+  CDGNO ( 0 : MAXPOG, 0: MAXPOG, KNDIMl +1) 

COMMON  / CRNTDM/  NMAX , MXNG , IQMAX , IQMAXl , IQMAX2 , IXCRNT , ICRNT , MXQEG , 

+  MXQOG 
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COMMON  / GCONST/  DKH , DKA , HH , HA , HSQ , DASQ , HSQN , DASQN , HHSQ , HDASQ , 

+  RAH , RAHSQ , DBA , DAH 

COMMON  / SUMLMT /  DHMAX , DAMAX , DSNG , IHMAX , IAMAX , ISNG 
COMMON  /NCONST/  DN , DNH , N 
COMMON  /GPQTMP/  CGNE , CDGNE , CGNO , CDGNO 
CHARACTER  FILEEVEN1*12 ,  FILEODDl*12 
SAVE  /GPQTMP/ , /GCONST/ , /SUMLMT/ 

C 

DKH = DKH IN 
DKA=DKAIN 

C  Initialize  the  matrix  CGNE,  CGNO,  CDGNE,  CDGNO 
DO  105  IC=  1 , KNDIM1+1 

DO  102  IB=  0 , MAXPEG 
DO  101  IA=  0 , MXPEG 
CGNE ( IA, IB , IC) =CZERO 
CDGNE {IA, IB, IC) =CZERO 

101  CONTINUE 

102  CONTINUE 

DO  104  ID=0 , MAXPOG 
DO  103  IE=0, MAXPOG 
CGNO (IE, ID, IC) =CZERO 
CDGNO ( IE , ID , IC ) =CZERO 

103  CONTINUE 

104  CONTINUE 

105  CONTINUE 
C 

C  Determine  the  maximum  number  of  terms  for  r-  and  m-  series  sums  and 
C  pass  them  through  /SUMLMT/: 

C 

REFC=FZERO*LOG (TWO) -LOG (PI2 ) -ONE 
REFH=HALF* (REFC-LOG (DKH) ) 

REFA=HALF* (REFC-LOG (DKA) ) 

C 

I QMX2  =  I QMAX+  2 
DHMAX= I QMX2 

DO  1100  WHILE  ( (DHMAX+HALF) * (LOG (DHMAX /DKH) +ONEN)  . LT .  REFH) 
DHMAX = DHMAX + ONE 
1100  CONTINUE 

IHMAX=INT (DHMAX) 

C 

DAMAX =AINT (DKA) +ONE 

DO  1200  WHILE  ( (DAMAX+HALF) * (LOG (DAMAX/DKA) +ONEN)  .LT.  REFA) 
DAMAX = DAMAX + ONE 
1200  CONTINUE 

C 

I AMAX= INT ( DAMAX ) 

C 

ISNG=INT (TWO*DKH+ONE-EPS8 ) 

ISNG=MAX ( IFPBIT , ISNG , IQMX2 ) 

C  Checking  dimensions. 

IF (IHMAX  .GT.  MXMREG)  THEN 

WRITE ( * , * )  ' Warning:  IHMAX  =  1  , IHMAX,  1  >  MXMREG  =  \  MXMREG 

WRITE (*,*)  ’IHMAX  IS  SET  TO  MXMREG  IN  XPQINI1 ’ 

IHMAX=MXMREG 
END  IF 

ISN1=MXMSNG-H-1 

IF (ISNG  .GT.  ISN1)  THEN 

WRITE ( * , * )  'Warning:  ISNG  =  1 , I SNG , 1  >  MXMSNG-N-1  =  ’ , ISN1 
WRITE (*,*)  'ISNG  IS  REDUCED  TO  MXMSNG-N-1  IN  XPQINIl ' 

ISNG=ISN1 
END  IF 
C 

DSNG=ISNG 

C 

FILEEVEN1= ' E 
FILEODDl=  * O 
NC1=INT (DKA) 

DC=NC1 
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NC=DKH*100000 . +NC1 
DC1=DKA-DC 
NC2=INT (1000 . *DC1) 

IGE=8*2* (MAXPEG+l) * (MAXPEG+2) 

IGO=8*2* (MAXPOG+1) * (MAXPOG+2) 

WRITE  (FILEEVEN1 (2:8) ,  1 (17.7) ')  NC 
WRITE  (FILEODD1 (2:8) ,*  (17.7) ')  NC 
WRITE  {FILEEVEN1 (10:12)  (13.3)  '  )  NC2 

WRITE  (FILEODD1 (10 : 12) ,  '  (13 . 3)  ' )  NC2 

OPEN  (28,ACCESS= 'DIRECT' , FILE=FILEEVEN1 , RECL=IGE , STATUS= ' UNKNOWN ' ) 
OPEN  (29, ACCESS= ' DIRECT ' , FILE = FI LEODDl , RECL= I GO, STATUS =' UNKNOWN ' ) 

C  Initialize  CGNE,  CDGNE ,  CGNO  and  CDGNO 
DO  1900  NI=1 , MXNG+1 

C  The  following  values  N,  DN  and  DNH  are  passed  to  the  G-computation 

C  related  subroutines  through  the  common  block  /NCONST/ . 

N=NI  -1 
DN=N 

DNH =DN+ HALF 
IRECE=0 

DO  1500  IQE=0 , MXQEG 
IQ=2*IQE 

DO  1300  IPE=0 , IQE-1 
IP=2*IPE 

CGNE ( IPE , IQE , NI ) =CGNE ( IQE ,  IPE ,  NI ) 

CDGNE ( IPE , IQE , NI ) = CDGNE ( IQE , IPE , NI ) 

1300  CONTINUE 

DO  1400  IPE=IQE, MXQEG 
IRECE=IRECE+1 
IP=2*IPE 

CALL  GDN (IP, IQ, GI , GDI , GR, GDR) 

GNE ( 1 , IRECE) =GR 

GNE ( 2 , IRECE ) =GI 

GNE ( 3 , IRECE) =GDR 

GNE ( 4 , IRECE) =GDI 

CGNE ( I PE , IQE , NI ) =DCMPLX ( GR ,  GI ) 

CDGNE ( I PE , IQE , NI ) =DCMPLX ( GDR , GDI ) 

1400  CONTINUE 

1500  CONTINUE 

WRITE  (28, REC=NI )  GNE 
IRECO=0 

DO  1800  IQO=0 , MXQOG 
IQ=2*IQO+l 

DO  1600  IPO=0,IQO-1 
IP=2  *IPO+l 

CGNO ( IPO , IQO , NI ) =CGNO ( IQO , IPO , NI ) 

CDGNO ( IPO , IQO , NI ) = CDGNO { IQO , IPO , NI ) 

1600  CONTINUE 

DO  1700  IPO=IQO, MXQOG 
IRECO=IRECO+l 
IP=2*IPO+l 

CALL  GDN ( IP, IQ, GI , GDI , GR, GDR) 

GNO ( 1 , IRECO) =GR 

GNO ( 2 , IRECO ) =GI 

GNO ( 3 , IRECO ) =GDR 

GNO (4, IRECO) =GDI 

CGNO ( I PO , IQO , NI ) =DCMPLX ( GR , GI ) 

CDGNO ( IPO , IQO , NI ) =DCMPLX ( GDR , GDI ) 

1700  CONTINUE 

1800  CONTINUE 

WRITE  (29 , REC=NI )  GNO 
1900  CONTINUE 
CLOSE (28) 

CLOSE (29) 

RETURN 

END 

C 

SUBROUTINE  GDN ( IPIN, IQIN, GIOUT , GDIOUT , GROUT , GDROUT) 

c********************************************************************** 
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anno  o  o 


C  This  subrouitne  sets  up  P  and  Q  dependent  constants  and  passes  along 
C  /PCONST/,  then  calls  subroutines  GDREG  and  GDSNG  to  compute  G(P,Q,N) 

C  and  its  derivative. 

Q* ********************************************************************* 

c 

INCLUDE  ' REALTP . INC  1 

COMMON  /PCONST/  DP , DQ, DS, DD, DSH, DDH, DSSQ, DDSQ , IP, IQ, IS,  ID 
C  Check  and  transform  input  variables. 

IP=IPIM 

IQ=IQIN 

ISPQ=IP+IQ 

IS=ISPQ/2 

IF  ( (ISPQ+1 ) /2  .GT.  IS)  THEN 
GIOUT=ZERO 
GDIOUT=ZERO 
GROUT=ZERO 
GDROUT = Z  E  RO 
WRITE  (*,  9001)  ISPQ 

9001  FORMAT ('The  parameters  P  and  Q  have  a  sum  of  the  ODD  integer  ',14, 
+/,',  G(P,  Q,  N)  and  its  derivative  have  been  set  to  0.') 

RETURN 
END  IF 

IF  (IP  .LT.  IQ)  THEN 
IP=IQIN 
IQ=IPIN 
END  IF 
C 

DP=IP 

DQ=IQ 

C 

ID=IS-IQ 

DS=IS 

DD=ID 

DSH=DS+HALF 
DDH = DD + HALF 
DSSQ=DS*DS 
DDSQ=DD*DD 
C 

CALL  GDREG (GI, GDI, GRR,GDRR) 

CALL  GDSNG ( GRS , GDRS ) 

GIOUT=GI 
GDIOUT=GDI 
GROUT=GRR+GRS 
GDROUT =GDRR+ GDRS 
RETURN 
END 
CC 

SUBROUTINE  GDREG (GI OUT, GDI OUT, GROUT, GDROUT) 

c********************************************************************** 

C  This  subrouitne  computes  the  regular  part  of  G(P,Q,N)  and  its 
C  derivative. 

q*  ************************************************  ********.**.**.****.***** 

C 

INCLUDE  1 REALTP . INC ' 

INCLUDE  1 GPQNDM . INC  1 

COMMON  / GCONST/  DKH, DKA, HH , HA, HSQ , DASQ , HSQN, DASQN, HHSQ, HDASQ , 

+  RAH ,  RAHSQ ,  DHA ,  DAH 

COMMON  /NCONST/  DN , DNH , N 

COMMON  /PCONST/  DP, DQ, DS, DD, DSH, DDH, DSSQ, DDSQ, IP, IQ, IS, ID 
COMMON  / SUMLMT /  DHMAX , DAMAX , DSNG , IHMAX , I AMAX , I SNG 
SAVE  /GCONST/, /SUMLMT/ 


Reserve  working  space  to  store  r-independent  numbers. 
DIMENSION  GIM(MXMREG) , GRRM (MXMREG) 

Computation  starts . 

Compute  overal  constant  factors: 


86 


GRRF=QUAR*DKH/PISQ/  (DSSQ-QUAR)  /  (QUAR-DDSQ)  /  (DN+ONE) 

GIF=HALF/DSH 

DJN=ZERO 

DO  300  JN=1 , N 
D  JN = D  JN  +  ONE 
HATN=HA/DJN 
GRRF=GRRF *  HATN*HATN 
GIF=GIF*HATN*HA/ (DJN+DSH) 

300  CONTINUE 

DJP=ZERO 

DO  400  JQ=1,IQ 
D  JP = D JP+ ONE 

GIF= (HHSQ/DJP/DJP) *GIF 
400  CONTINUE 

DO  500  JP=IQ+1 , IP 
DJP=DJP+ONE 
GIF= (HH/DJP) *GIF 
500  CONTINUE 

IF  { (ID+1) /2  .GT.  ID/2)  GIF--GIF 
C 

C  Compute  GI,  GDI ,  GRREG,  and  GDRREG . 

C  Compute  r-independent  factors  and  store  in  GIM  and  GRRM: 

SMH=DSH+ONEN 
SM1=DS 
SM2=DS+DS 
PM1=DP 
QM1=DQ 
THRH=HALF 
DMl=ZERO 
THRS=HALF+DS 
THRD=HAIiF+DD 
THRDN=HALF~DD 
THRSN=HALF-DS 

DO  600  JM=1 , IHMAX 
SMH=SMH+ONE 
SMl=SMl+ONE 
SM2  =  SM2  +  ONE 
PM1  =  PM1  +  ONE 
QM1=QM1+0NE 
THRH=THRH+ONE 
DMl=DMl+ONE 
THRS=THRS+ONE 
THRD=THRD+ONE 
THRDN = THRDN  +  ONE 
THRSN=THRSN+ONE 

GIM(JM)  =  (SMH/SM2)  *  (SMH/PMl)  *  (SM1/QM1)  *  (HSQN/DM1 ) 

GRRM  ( JM)  =  (THRH/THRS)  *  (HSQN/THRD)  *  (DM1 /THRDN)  *  (DM1/THRSN) 
600  CONTINUE 

C 

C  Compute  r-  and  m-sum  for  GI  and  GRR. 

C  Setup  initial  r  related  values 

DJR=DAMAX 
DNR=DN+DJR 
DNHR=DNR-HALF 
D2NR=DN+DNR 
DNSHR=DNR+DSH 
DN2R=DNR+ONE 

C  m-sum  for  r=IAMAX 

DN  SHRM= DNSHR+ DHMAX 
DN2  RM=  DN2  R+DHMAX 
ANNEXTR=ONE+GIM ( IHMAX) /DNSHRM 
ANR0R=ONE+GRRM( IHMAX) /DN2RM 
DO  700  JM= IHMAX- 1 , 1 , -1 
DNSHRM=DNSHRM+ONEN 
DN2RM=DN2RM+ONEN 

ANNEXTR=ONE+ANNEXTR*GIM ( JM) /DNSHRM 
ANR0R=ONE+ANR0R*GRRM(JM) /DN2RM 
700  CONTINUE 
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on  no  non  on 


ANNEXT=ANNEXTR 
DANNEXT = DNR  * ANNEXTR 
ANRO=ANROR 
DANR  0 =DNR  *  ANR  0  R 

DO  900  JR=IAMAX, 1,-1 
C  Compute  factors  for  r=JR 

ATMP= (DNHR/D2NR) * (DASQN/DJR) 

ATMP I -ATMP / DNSHR 
ATMPR= ATMP / DN2 R 

C  Setup  new  r  related  values  for  r=JR-l 

DJR=DJR+ONEN 
DNR= DNR + ONEN 
DNHR=DNHR+ ONEN 
D2NR=D2NR+ONEN 
DNSHR=DNSHR+ONEN 
DN  2  R = DN2  R + ONEN 
C  m-sum  for  r=JR-l 

DNSHRM=DNSHR+DHMAX 
DN2 RM= DN2 R+DHMAX 
ANNEXTR=ONE+GIM ( IHMAX) /DNSHRM 
ANROR=ONE+GRRM ( IHMAX) /DN2RM 
DO  800  JM= IHMAX- 1 ,1,-1 
DNSHRM=DNSHRM+ONEN 
DN2 RM= DN2 RM+ ONEN 

ANNEXTR=ONE+ANNEXTR*GIM(JM) /DNSHRM 
ANROR=ONE+ANROR*GRRM ( JM) /DN2RM 
800  CONTINUE 

ANNEXT=ANNEXTR+ANNEXT * ATMPI 
DANNEXT = DNR  *  ANNEXT  R + DANNEXT  *  ATMP  I 
ANR0=ANR0R+ANR0*ATMPR 
DANR  0 = DNR  *  ANR  0  R + DANR  0  *  ATMPR 
900  CONTINUE 

C 

GIOUT=ANNEXT*GIF 
GDIOUT=DANNEXT*GIF 
GROUT= ANR  0  *  GRRF 
GDROUT=DANR0  *GRRF 
C 

RETURN 

END 

C 

SUBROUTINE  GDSNG (GROUT , GDROUT) 

c********************************************************************** 

C  This  subrouitne  computes  the  'singular'  part  of  G(P,Q,N)  and  its 
C  derivative. 

£*******************************************★************************** 

C 

INCLUDE  1 REALTP . INC ' 

INCLUDE  ' GPQNDM . INC ' 

COMMON  / GCONST/  DKH , DKA , HH , HA , HSQ , DASQ , HSQN , DASQN , HHSQ , HDASQ , 

+  RAH, RAHSQ, DHA, DAH 

COMMON  /NCONST/  DN , DNH , N 

COMMON  /PCONST/  DP, DQ, DS, DD, DSH, DDH, DSSQ,DDSQ, IP, IQ,  IS,  ID 
COMMON  /SUMLMT/  DHMAX, DAMAX, DSNG, IHMAX, IAMAX, ISNG 
SAVE  /GCONST/, /SUMLMT/ 


Reserve  working  space  to  store  r-independent  numbers. 
DIMENSION  GS1M (MXMSNG) , GS2M (MXMSNG) , 

+  GR2M (MXMSNG) , GR2R ( 0 : MXMSNG) , GDR2 R ( 0 : MXMSNG ) 

Computation  starts. 

MXGR2M=ISNG+N 

Compute  overal  constant  factors: 

GRSF=QUAR/PISQ/DKH 

Compute  GR1 ,  GDR1 ,  GR2 ,  and  GDR2 : 
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c 


Compute  initial  constants: 


1000 


1100 


1200 


C 

c 


1300 

C 

C 


S1LHA=TW0*L0G (HXD/RAH) 

DK=HALF 

S1SDO=ZERO 

DO  1000  JK=1 , ID 
S1SDO=S1SDO+ONE/DK 
DK=DK+ONE 
CONTINUE 
S1SDO=TWO*S1SDO 

DO  1100  JK=ID+1,  IS 
SlSD0=SlSD0+ONE/DK 
DK=DK+ONE 
CONTINUE 

S1SDO=-TWO*S1SDO 

SN0N=ZERO 

SNIO^ZERO 

SN20=ZERO 

IF  (N  .GT.  0)  THEN 
DJN=ONE 
DJNH=HALF 

DO  1200  JN=1 , N-l 
D JNI = ONE / D JN 
D JNHI =ONE / D JNH 
SN1 0 = D JNI *  D JNHI  +  SN1 0 

SN2 0 =D JNHI *D JNHI*  (  (ONE-QUAR*DJNI )  * HALF * D JNI + ONE)  +SN20 
DJN=DJN+ONE 
D JNH = D JNH + ONE 
CONTINUE 
D JNI = ONE / D  JN 
D  JNH  I = ONE  /  D  JNH 

SN0N=S1LHA+HALF*DJNHI-SN10*QUAR 
SN10= { D  JNI *D JNHI  +  SN1 0 ) *QUAR 

SN20=HALF*  (DJNHI*DJNHI*  (  (QUAR*DJNI+ONEN)  * HALF *D  JNI +ONEN)  -SN20) 
END  IF 

SN10=S1LHA-SN10 

SN20=PISQ*THIR+SN20 


Compute  r-independent  terms  and  store  in  GR2M,  GS1M,  and  GS2M: 

GR2 JM=ONE 

S1SDM=S1SD0 

S2SDM=ZERO 

DM=ZERO 

DMH=-HALF 

DO  1300  JM=1 , MXGR2M 
DM=DM+ONE 
DMH = DMH + ONE 
DMI=ONE/DM 
DMH I = ONE / DMH 
DMHSQ=DMH*DMH 
DMHSS=DMHSQ-DSSQ 
DMHSD=DMHSQ~DDSQ 

GR2  JM=  (DMHSS*DMI*DMI )  *  (DMHSD*DMI*DMHI )  *GR2  JM 
GR2M { JM) =GR2JM 
TWSMHS=TWO*DSSQ/DMHSS 
TWDMHD=TWO  *DDSQ / DMHSD 

S1SDM=S1SDM- (TWSMHS+TWDMHD+DMI+ONE)  *DMHI 
S2SDM=S2SDM-  (  (TWSMHS+THR)  *TWSMHS+  (TWDMHD+THR)  *TWDMHD+ 

+  (ONE-QUAR*DMI) *TWO*DMI+ONE) /DMHSQ 

GS1M ( JM) =S1SDM 
GS2M ( JM) =S2SDM 
CONTINUE 

Compute  m-sum  for  GRl ,  GDR1  and  store  in  GR2R,  GDR2R 
IF  (N  .GT.  0)  THEN 
D JR = DN +  ONEN 
DJRH=DJR-HALF 
SN0R=SN0N 
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1400 


1500 

1600 

C 

c 


1700 


SGR1M=SN0R+S1SD0 
SDGR1M=D JR*  SGR1M+ONEN 
DJM=ONE 
DJRM=DJR 

DO  1400  JM-1 , N-l 
SR1M=SN0R+GS1M ( JM) 

HSQF=HSQ/DJM/DJRM 
SGR1M=SGR1M*HSQF+SR1M*GR2M ( JM) 
SDGR1M=SDGR1M*HSQF+ (DJR*SR1M+0NEN) *GR2M(JM) 
D JM = D JM + ONE 
D JRM= D JRM+ OMEN 
CONTINUE 

GR2R {N-l ) =SGRlM*TWO 
GDR2R {N-l ) =SDGRlM*TWO 

DO  1600  JR=N-2 ,0,-1 

SN  0  R= SN 0  R+ ONE / D JRH + ONE / (DN-DJR) -ONE/ (DN+DJR) 

DJR=DJR+ONEN 

DJRH=DJRH+ONEN 

SGR1M=SN0R+S1SD0 

SDGR1M=DJR*  SGR1M+ONEN 

DJM=ONE 

DJRM=DJR 

DO  1500  JM=1 , JR 
SR1M=SN0R+GS1M { JM) 

HSQF=HSQ/DJM/DJRM 
SGRlM=SGRlM*HSQF+SRlM*GR2M( JM) 
SDGRlM=SDGRlM*HSQF+ (DJR*SRlM+ONEN) *GR2M(JM) 

D JM = D JM+  ONE 
DJRM=DJRM+ONEN 

CONTINUE 

GR2R (JR) =SGRlM*TWO 
GDR2R( JR) =SDGRlM*TWO 
CONTINUE 
END  IF 


Compute  m-sum  for  GR2 ,  GDR2  and  store  in  GR2R,  GDR2R 
DRN=DM 
SNlR=SN10 
SN2R=SN20 
SR2MM=SN1R+S1SD0 
SR2RMS=SR2MM*SR2MM+SN2R 
DR2RMS=DRN*SR2 RMS-TWO* SR2MM 
DJM=ONE 
DJRM^DRN 

DO  1700  JM= 1 , N 
SR2MM=SN1R+GS1M ( JM) 

SR2RM=SR2MM*SR2MM+SN2R+GS2M ( JM) 
DR2RM=DRN*SR2RM-TWO*SR2MM 
HSQF=HSQ/DJM/DJRM 
SR2RMS=SR2RMS*HSQF+SR2RM*GR2M ( JM) 
DR2RMS=DR2RMS*HSQF+DR2RM*GR2M ( JM) 

D JM = D JM+ ONE 
DJRM= DJRM+ ONEN 
CONTINUE 
GR2R (N) =SR2RMS 
GDR2R (N) =DR2RMS 
DJR=ZERO 

DO  1900  JR=N+1 , MXGR2M 
D JR = D JR  +  ONE 
DRN = DRN +  ONE 
D JRI = ONE / D JR 
DRNHI=ONE/ (DRN- HALF) 

DR2NI=ONE/ (DRN+DN) 

DNHRF=DNH*DRNHI*DR2NI 

SN1R=DJRI-DNHRF+SN1R 

SN2R=DJRI*DJRI- (DRNHI+DR2NI ) *DNHRF+SN2R 

SR2MM=SN1R+S1SD0 

SR2RMS=SR2MM* SR2MM+SN2R 
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1800 


1900 

c 

c 

c 

c 

c 

c 


2000 


2100 

C 

c 


2200 

2300 


DR2RMS=DRN*SR2RMS-TWO*SR2MM 

DJM=ONE 

DJRM=DRN 

DO  1800  JM=1 , JR 
SR2MM=SN1R+GS1M ( JM) 

SR2RM=SR2MM*SR2MM+SN2R+GS2M( JM) 

DR2RM=DRN*SR2RM-TWO*SR2MM 
HSQF=HSQ / D JM / DJRM 
SR2RMS=SR2RMS*HSQF+SR2RM*GR2M ( JM) 
DR2RMS=DR2RMS*HSQF+DR2RM*GR2M ( JM) 

D JM= D JM+ ONE 
DJRM=DJRM+ONEN 
CONTINUE 
GR2R (JR) =SR2RMS 
GDR2R ( JR) =DR2RMS 
CONTINUE 

Compute  r-sum: 

IF  (RAHSQ  .GT.  HALF)  THEN 

Convergence  acceleration  via  Euler-Abel  Transfomation 

DJR=ZERO 
DJRH=DJR-HALF 
DJRN=DJR+DN 
D JRNN =  D JR - DN 

IF  {N  .GT.  0)  THEN 
RFCTOR=ONE / DN 
GR2R ( 0 ) =RFCTOR*GR2R ( 0 ) 

GDR2R ( 0 ) =RFCTOR*GDR2R ( 0 ) 

DO  2000  JR=1 , N-l 
DJR=DJR+ONE 
DJRH=DJRH+ONE 
DJRN=DJRN+ONE 
D  JRNN =D  JRNN + ONE 

RFCTOR= ( D JR *  D JRH / D JRN / D JRNN ) *RAHSQ*RFCTOR 
GR2R { JR) =RFCTOR*GR2R ( JR) 

GDR2R (JR) =RFCTOR*GDR2R ( JR) 

CONTINUE 
DJR=DJR+ONE 
DJRH=DJRH+ONE 
D  JRN = D  JRN + ONE 
DJRNN= DJRNN+ ONE 

RFCTOR=HALF* (HALF-DN) *RAHSQ*RFCTOR 
GR2R (N) =RFCTOR*GR2R (N) 

GDR2R (N) =RFCTOR*GDR2R (N) 

ELSE 

RFCTOR=ONE 

END  IF 

DO  2100  JR=N+1 , MXGR2M 
DJR=DJR+ONE 
D JRH = D JRH + ONE 
D  JRN = D  JRN + ONE 
DJRNN=  DJRNN+ ONE 

RFCTOR= ( D JR *  D JRH / D JRN / D JRNN ) * RAHSQ* RFCTOR 
GR2R (JR) =RFCTOR*GR2R (JR) 

GDR2R ( JR) =RFCTOR*GDR2R (JR) 

CONTINUE 

Compute  zeroth-order  delta-r  term: 

DO  2300  JDELTR=1,MXGR2M 

DO  2200  JR=MXGR2M, JDELTR, -1 
GR2R ( JR) =GR2R { JR-1 ) -GR2R ( JR) 

GDR2R ( JR) =GDR2R ( JR-1 ) -GDR2R (JR) 

CONTINUE 

CONTINUE 

SGR2R=HALF*GR2R (MXGR2M) 
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ooh-^n  to  to  non 


SDGR2  R=HALF  *GDR2  R (MXGR2M) 

DO  2400  JR=MXGR2M-1 ,0,-1 
SGR2R= ( SGR2R+GR2R (JR) ) *HALF 
SDGR2R= ( SDGR2R+GDR2R ( JR) ) *HALF 
2400  CONTINUE 

GROUT=GRSF*SGR2R 
GDROUT=GRSF*  SDGR2 R 
C 

ELSE 

Direct  summation  when  RAHSQ  is  no  greater  than  1/2: 

DJR-DSNG 
DJRN=DJR+DN 
DJRNH=DJRN-HALF 
DJR2N=DJRN+DN 

RFCTOR= ( D JRN*  D JRNH / D JR / D JR2N ) * RAHSQ 
SGR2R=RFCTOR*GR2R (MXGR2M) 

SDGR2R=RFCTOR*GDR2R (MXGR2M) 

DO  2500  JR=MXGR2M-1,N+1, -1 
DJR=DJR+ONEN 
DJRN=DJRN+ONEN 
DJRNH=DJRNH+ONEN 
DJR2N=DJR2N+ONEN 

RFCTOR= ( D JRN *  D JRNH / D JR / D JR2 N ) * RAHSQ 
SGR2R= (GR2R( JR) -SGR2R) *RFCTOR 
SDGR2R= (GDR2R{ JR) -SDGR2R) *RFCTOR 
00  CONTINUE 

IF  (N  .EQ.  0)  THEN 
GROUT = (GR2R (0) -SGR2R) *GRSF 
GDROUT= (GDR2R ( 0 ) -SDGR2R) *GRSF 
ELSE 

RFCTOR=HALF* (HALF-DN) * RAHSQ 
SGR2R= (GR2R (N) -SGR2R) *RFCTOR 
SDGR2R= (GDR2R (N) -SDGR2R) *RFCTOR 
DJR=DN 

D JRH = D JR - HAL F 
DJRN=DJR+DN 
D JRNN = D JR - DN 

DO  2600  JR=N-1 ,1,-1 
DJR=DJR+ONEN 
D JRH = D JRH + ONEN 
D JRN= D JRN+  ONEN 
D JRNN= D JRNN+ ONEN 

RFCTOR=  ( D  JR  *  D  JRH  /  D  JRN  /  D  JRNN )  *  RAHSQ 
SGR2R= (GR2R (JR) -SGR2R) *RFCTOR 
SDGR2R= (GDR2R (JR) -SDGR2R) *RFCTOR 
00  CONTINUE 

GROUT= (GR2R ( 0 ) -SGR2R) *GRSF/DN 
GDROUT= (GDR2R ( 0 ) -SDGR2R) *GRSF/DN 
END  IF 
END  IF 
RETURN 
END 


SUBROUTINE  BCKSCFL 

SUBROUTINE  BCKSCFL ( IR , CESTH , CESPH ) 

★  ★★★★it********************************************************************** 

INCLUDE  ' REALTP . INC ' 

INCLUDE  ‘ CMPXTP . INC  1 
INCLUDE  ' MAINDM . INC ' 

C 

REAL* 4  AKPHPR , AKPHPI , AKZPR , AKZPI , ALPHPR , ALPHPI , ALZPR , ALZPI , 
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+  AKPHNR , AKPHNI , AKZNR , AKZNI , ALPHNR , ALPHNI , ALZNR , ALZNI 

DIMENSION  CIW (KCRNT) , CKLN (KCRNT) 

DIMENSION  CIWO ( KXCRT ) , CKLNO ( KXCRT ) , CXPQNO ( KXCRT , KXCRT ) 

DIMENSION  CRPQ1 (KCRNT, KCRNT) , CRPQ2 ( KCRNT , KCRNT ) , 

+  CXPQN ( KCRNT , KCRNT ) , CXRPQ ( KCRNT , KCRNT ) 

DIMENSION  CKPHP(61, 61) , CKZP <61, 61) , CLPHP ( 61 , 61 )  , 

+  CLZP (61,61) ,  CKPHN (61,61), CKZN (61,61)  , 

+  CLPHN (61,61) , CLZN ( 61 , 61 ) 

COMMON  / INPUT2 /  IE, IM, THETAI , THESINI , THECOSI , RTHEI 

COMMON  / INPUT4 /  IZ , IK , IS , NYSM 

COMMON  /NCONST/  DN , DNH , N 

COMMON  /XPQTMP/  CXPQN, CXRPQ , CXPQNO 

COMMON  / XPQTMP 1 /  CRPQ1,CRPQ2 

COMMON  /CKLMTX/  CKPHP, CKZP, CLPHP , CLZP, CKPHN, CKZN, CLPHN, CLZN 
COMMON  / CRNTDM/  NMAX, MXNG, IQMAX, IQMAX1 , IQMAX2 , IXCRNT, ICRNT, MXQEG, 

+  MXQOG 

COMMON  / INPUT 3 /  RTHE , RDELT , RPHI , RDELP , NTHTAO , NPHI , THESIN , THECOS , 

+  RHPI 

CESTH=CZERO 
CESPH=CZERO 
DO  100  JQ=1, KXCRT 
CKLNO (JQ)=CZERO 
100  CONTINUE 

DO  200  JQ=1, KCRNT 
CKLN ( JQ) =CZERO 
200  CONTINUE 

DO  400  IF=-30 , 3 0 
DO  300  JZ=-3 0,30 
CKPHP ( IF , JZ ) =CZERO 
CKZP (IF, JZ)=CZERO 
CLPHP (IF, JZ)=CZERO 
CLZP (IF , JZ) =CZERO 
CKPHN ( IF , JZ ) =CZERO 
CKZN ( IF , JZ ) =CZERO 
CLPHN (IF, JZ)=CZERO 
CLZN ( IF , JZ ) =CZERO 
3  00  CONTINUE 

400  CONTINUE 
C 

IF  (RTHEI  .EQ.  0 . DO )  GO  TO  2000 
C  If  incident  angle  is  not  equal  to  0,  use  this  loop 
DO  1800  NA=0 , NMAX 
C 

CALL  INCIDNT(NA, CIWO, CIW) 

C 

IF  (NA  .EQ.  0)  THEN 
CALL  XPQ0 
ELSE 

CALL  XPQN(NA) 

END  IF 

C  If  the  cylinder  is  made  of  perfect  conductor  and  no  coating  on  it 
IF  (IZ  .EQ.  0)  THEN 
DN=NA 

C  Use  IMSL  library  to  solve  linear  system 

CALL  DLSACG ( IXCRNT , CXPQNO , KXCRT , CIWO , 1 , CKLNO ) 

CALL  ESCFAR (NA, CKLNO , CKLN, CETHN, CEPHN) 

C 

C  If  we  calculate  in  X  and  Y  components  as  theta  approaches  0  or  pi, 

C  then  theta  component  is  cosin  phi  in  X  direction  plus  sine  phi  in  Y 
C  direction,  and  phi  component  is  negative  sine  phi  in  X  direction  plus 
C  cosin  phi  in  Y  direction. 

IF  (RTHE  .EQ.  0.D0)  THEN 
CETHN=COS (RPHI ) *CETHN+SIN (RPHI) *CEPHN 
CEPHN=-SIN (RPHI ) *CETHN+COS (RPHI ) *CEPHN 
ELSEIF  (RTHE  . EQ .  PI)  THEN 
CETHN=-COS (RPHI ) *CETHN-SIN (RPHI ) *CEPHN 
CEPHN=-SIN (RPHI ) *CETHN+COS (RPHI ) *CEPHN 
END  IF 
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c 

CESTH=CESTH+CETHN 

CESPH=CESPH+CEPHN 

IF  (NA  .NE.  0)  THEN 
IF  (IE  .EQ.  1)  THEN 
DO  600  1=2 , KXCRT , 2 
CKLNO (I) =-CKLN0 (I) 

600  CONTINUE 

END  IF 

IF  (IM  .EQ.  1)  THEN 
DO  800  1=1, KXCRT, 2 
CKLNO (I) = -CKLNO (I) 

800  CONTINUE 

END  IF 

NA1=-NA 

DN=NA1 

CALL  ESCFAR (NA1 , CKLNO , CKLN, CETHN, CEPHN) 

IF  (RTHE  .EQ.  0 . DO )  THEN 
CETHN=COS (RPHI ) *CETHN+SIN (RPHI ) *CEPHN 
CEPHN=-SIN (RPHI ) *CETHN+COS (RPHI) *CEPHN 
ELSEIF  (RTHE  .EQ.  PI)  THEN 
CETHN=-COS (RPHI) *CETHN-SIN (RPHI) * CEPHN 
CEPHN=-SIN (RPHI ) *CETHN+COS (RPHI) * CEPHN 
END  IF 
C 

CESTH=CESTH+CETHN 
CESPH=CESPH+CEPHN 
END  IF 
END  IF 

C  If  the  cylinder  is  coated  with  anisotropic  material 
IF  (IZ  .EQ.  1)  THEN 
DO  1100  1=1 , KCRNT 
DO  1000  J=l, KCRNT 
CXRPQ ( I , J) =CXPQN ( I , J) +CRPQ1 (I , J) 

1000  CONTINUE 
1100  CONTINUE 
DN=NA 

CALL  DLSACG ( ICRNT , CXRPQ , KCRNT , CIW , 1 , CKLN ) 

IF  ((IK  .EQ.  1)  .AND.  (IR  . EQ .  45))  THEN 
CALL  KLCRNT(DN, CKLN, CIW) 

END  IF 

CALL  ESCFAR (NA, CKLNO , CKLN, CETHN, CEPHN) 

IF  (RTHE  .EQ.  0.D0)  THEN 
CETHN=COS (RPHI ) *CETHN+SIN (RPHI ) *CEPHN 
CEPHN=-SIN (RPHI ) *CETHN+COS ( RPHI ) *CEPHN 
ELSEIF  (RTHE  .EQ.  PI)  THEN 
CETHN=-COS  (RPHI)  *  CETHN- SIN  (RPHI)  * CEPHN 
CEPHN= - SIN ( RPHI ) *  CETHN+COS ( RPHI ) *  CEPHN 
END  IF 
C 

CESTH=CESTH+ CETHN 
CESPH=CESPH+CEPHN 

IF  (NA  .NE.  0)  THEN 

IF  (NYSM  .EQ.  0)  THEN 
DO  1300  1=1, KCRNT 
DO  1200  J=l, KCRNT 
CXRPQ ( I , J) =CXPQN ( I , J) +CRPQ2 (I , J) 

1200  CONTINUE 

1300  CONTINUE 

ELSE 

DO  1500  1=1, KCRNT 
DO  1400  J=l, KCRNT 
CXRPQ ( I , J) =CXPQN ( I , J) +CRPQ1 ( I , J) 

1400  CONTINUE 

1500  CONTINUE 

END  IF 

CALL  DLSACG ( ICRNT , CXRPQ , KCRNT, CIW, 1 , CKLN) 

IF  (IE  .EQ.  1)  THEN 
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DO  1600  I=2,KCRNT,4 
11=1+1 

CKLN ( I ) = -CKLN { I ) 

CKLN ( I 1 ) = -CKLN (II) 

1600  CONTINUE 

END  IF 

IF  (IM  .EQ.  1)  THEN 
DO  1700  1=1 , KCRNT , 4 
12=1+3 

CKLN ( I ) =  --CKLN ( I ) 

CKLN (12 ) =-CKLN ( 12 ) 

1700  CONTINUE 

END  IF 

NAl=-NA 

DN=NA1 

IF  ((IK  .EQ.  1)  .AND.  (IS  . EQ .  199))  THEN 
CALL  KLCRNT (DN, CKLN, Cl W) 

END  IF 

CALL  ESCFAR (NAl , CKLNO , CKLN, CETHN, CEPHN) 

IF  (RTHE  .EQ.  0.D0)  THEN 
CETHN=COS (RPHI ) *CETHN+SIN (RPHI ) *CEPHN 
CEPHN=- SIN (RPHI ) *CETHN+COS (RPHI ) *CEPHN 
ELSEIF  (RTHE  .EQ.  PI)  THEN 
CETHN=-COS (RPHI ) *CETHN-SIN (RPHI ) *CEPHN 
CEPHN=-SIN (RPHI ) *CETHN+COS (RPHI ) *CEPHN 


CESTH=CESTH+CETHN 
CESPH=CESPH+CEPHN 
END  IF 
END  IF 

1800  CONTINUE 

C  If  the  cylinder  is  coated  with  anisotropic  material  and  we  want  to  calculate 
C  equivalent  currents  K  nad  L  on  inner  and  outer  surfaces. 

IF  (IK  .EQ.  1)  THEN 

IF  (IS  .EQ.  199)  THEN 
OPEN (31, FILE= ' khzz41 . dz ' , STATUS = ‘ UNKNOWN ' ) 

OPEN (32 , FILE= ’ lhzz41 . dz 1 , STATUS =* UNKNOWN ' ) 

END  IF 

DO  1900  IF=-30 , 30 
DO  1850  JZ=-30, 30 
AKPH  PR = REAL ( CKPHP ( IF , JZ ) ) 

AKPHPI=IMAG ( CKPHP ( IF , JZ ) ) 

AKZPR=REAL (CKZP ( IF, JZ ) ) 

AKZPI=IMAG (CKZP (IF, JZ) ) 

AK  PHNR = REAL ( CKPHN ( IF , JZ ) ) 

AKPHNI=IMAG (CKPHN (IF, JZ) ) 

AKZNR=REAL ( CKZN ( IF , JZ ) ) 

AKZNI=IMAG (CKZN (IF, JZ) ) 

ALPHPR=REAL ( CLPHP ( IF , JZ ) ) 

ALPHPI  =  IMAG (CLPHP (IF  ,  JZ ) ) 

ALZPR=REAL ( CLZP ( IF , JZ ) ) 

ALZPI=IMAG (CLZP ( IF, JZ) ) 

AL  PHNR = REAL ( CLPHN ( IF , JZ ) ) 

ALPHNI=IMAG ( CLPHN ( IF , JZ ) ) 

ALZNR=REAL ( CLZN ( IF , JZ ) ) 

ALZNI=IMAG (CLZN ( IF , JZ ) ) 

WRITE (31,*)  AKPHPR, '  ‘ ,AKPHPI 

WRITE (31,*)  AKZPR, 1  ‘ ,AKZPI 

WRITE ( 32 , * )  ALPHPR, 1  ’,ALPHPI 
WRITE (3 2,*)  ALZPR, '  ',ALZPI 
WRITE ( 31 , * )  AKPHNR, '  ' , AKPHNI 

WRITE (31,*)  AKZNR, '  * ,AKZNI 

WRITE (32 , * )  ALPHNR, '  ’ ,ALPHNI 

WRITE (3 2,*)  ALZNR, 1  ',ALZNI 
1850  CONTINUE 

1900  CONTINUE 
CLOSE (31) 
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CLOSE (32) 

END  IF 
RETURN 

C  If  incident  angle  is  equal  to  zero,  the  program  should  go  to  this  loop 
C  because  it  need  to  compute  when  n=+l  and  n=~l  only. 

2000  NA=1 
DN=NA 

CALL  INCIDNT (NA, CIWO , CIW) 

CALL  XPQN(NA) 

C 

IF  (12  .EQ.  0)  THEN 

CALL  DLSACG ( IXCRNT , CXPQNO , KXCRT , CIWO , 1 ,  CKLNO ) 

CALL  ESCFAR (NA, CKLNO , CKLN, CETHN, CEPHN) 

C 

C  If  we  calculate  in  X  and  Y  components  as  theta  approaches  0  or  pi, 

C  then  theta  component  is  cosin  phi  in  X  direction  plus  sine  phi  in  Y 
C  direction,  and  phi  component  is  negative  sine  phi  in  X  direction  plus 
C  cosin  phi  in  Y  direction. 

IF  (RTHE  .EQ.  0.D0)  THEN 
CETHN=COS (RPHI ) *CETHN+SIN (RPHI ) * CEPHN 
CEPHN=-SIN (RPHI ) *CETHN+COS (RPHI ) *CEPHN 
ELSEIF  (RTHE  .EQ.  PI)  THEN 
CETHN=-COS (RPHI) *CETHN-SIN< RPHI) * CEPHN 
CEPHN=-SIN (RPHI ) *CETHN+COS( RPHI) *CEPHN 
END  IF 
C 

CESTH=CESTH+CETHN 

CESPH=CESPH+CEPHN 

IF  (NA  .NE.  0)  THEN 

IF  (IE  .EQ.  1)  THEN 
DO  2400  1=2, KXCRT, 2 
CKLNO (I)=-CKLN0 (I) 

2400  CONTINUE 

END  IF 

IF  (IM  .EQ.  1)  THEN 
DO  2500  1=1, KXCRT, 2 
CKLNO (I) =-CKLN0 (I) 

2500  CONTINUE 

END  IF 
NA1=-NA 
DN=NA1 

CALL  ESCFAR (NAl , CKLNO , CKLN, CETHN, CEPHN) 

IF  (RTHE  .EQ.  0.D0)  THEN 
CETHN=COS (RPHI ) *CETHN+SIN (RPHI ) * CEPHN 
CEPHN=-SIN (RPHI ) *CETHN+COS (RPHI) *CEPHN 
ELSEIF  (RTHE  . EQ .  PI)  THEN 
CETHN=-COS (RPHI ) * CETHN- SIN (RPHI ) * CEPHN 
CEPHN=-SIN (RPHI ) *CETHN+COS (RPHI) * CEPHN 
END  IF 
C 

CESTH=CESTH+CETHN 
CESPH=CESPH+CEPHN 
END  IF 
END  IF 
C 

IF  (IZ  .EQ.  1)  THEN 
DO  3100  1=1 , KCRNT 
DO  3000  J=l, KCRNT 

CXRPQ ( I , J) =CXPQN ( I , J) +CRPQ1 ( I , J) 

3000  CONTINUE 
3100  CONTINUE 
DN=NA 

CALL  DLSACG ( ICRNT , CXRPQ , KCRNT , CIW , 1 , CKLN ) 

C 

IF  ((IK  .EQ.  1)  .AND.  (IR  . EQ .  45))  THEN 
CALL  KLCRNT(DN, CKLN, CIW) 

END  IF 


C 
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c 


3200 

3300 


3400 

3500 


3600 


3700 


C 


C 


C 


C 


CALL  ESCFAR (NA, CKLNO , CKLN, CETHN, CEPHN) 

IF  {RTHE  .EQ.  O.DO)  THEN 
CETHN=COS (RPHI) * CETHN+S IN { RPHI )  * CEPHN 
CEPHN=-SIN(RPHI) *CETHN+COS (RPHI)*CEPHN 
ELSEIF  (RTHE  . EQ .  PI)  THEN 
CETHN= -COS (RPHI) * CETHN- SIN (RPHI )* CEPHN 
CEPHN=-SIN (RPHI ) *CETHN+COS (RPHI) *CEPHN 
END  IF 

CESTH=CESTH+CETHN 

CESPH=CESPH+CEPHN 

IF  (NA  .NE.  0)  THEN 

IF  (NYSM  .EQ.  0)  THEN 
DO  3300  1=1 , KCRNT 
DO  3200  J=l, KCRNT 
CXRPQ ( I ,  J) =CXPQN (I , J) +CRPQ2 ( I , J) 

CONTINUE 

CONTINUE 

ELSE 

DO  3500  1=1 , KCRNT 
DO  3400  J=1 , KCRNT 
CXRPQ (I , J) =CXPQN (I, J) +CRPQ1 ( I , J ) 

CONTINUE 

CONTINUE 

END  IF 

CALL  DLSACG ( ICRNT , CXRPQ , KCRNT , CIW, 1 , CKLN) 

IF  (IE  .EQ.  1)  THEN 
DO  3600  1=2, KCRNT, 4 
11=1+1 

CKLN ( I ) = -CKLN ( I ) 

CKLN (II) = -CKLN (II) 

CONTINUE 

END  IF 

IF  (IM  .EQ.  1)  THEN 
DO  3700  1=1, KCRNT, 4 
12=1+3 

CKLN { I ) = -CKLN ( I ) 

CKLN (12 ) =-CKLN (12 ) 

CONTINUE 

END  IF 

NA1=-NA 

DN=NA1 

IF  ((IK  .EQ.  1)  .AND.  (IS  . EQ .  199))  THEN 
CALL  KLCRNT(DN, CKLN, CIW) 

END  IF 

CALL  ESCFAR (NA1 , CKLNO , CKLN , CETHN , CEPHN) 

IF  (RTHE  .EQ.  O.DO)  THEN 
CETHN=COS (RPHI) *CETHN+SIN (RPHI) *CEPHN 
CEPHN=-SIN (RPHI ) *CETHN+COS (RPHI) *CEPHN 
ELSEIF  (RTHE  .EQ.  PI)  THEN 
CETHN=-COS (RPHI) *CETHN-SIN (RPHI) *CEPHN 
CEPHN=-SIN (RPHI ) *CETHN+COS (RPHI) *CEPHN 
END  IF 

CESTH=CESTH+CETHN 
CESPH=CESPH+CEPHN 
END  IF 
END  IF 

IF  (IK  .EQ.  1)  THEN 

IF  (IS  .EQ.  199)  THEN 
OPEN (31, FILE= ' khzz41.dz 1 , STATUS =' UNKNOWN ' ) 
OPEN (32, FILE= 1 lhzz41.dz* , STATUS =' UNKNOWN ' ) 
END  IF 

DO  4200  IF=-30 , 30 
DO  4100  JZ=-30, 30 
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AKPHPR=REAL ( CKPHP ( IF , JZ ) ) 

AKPHPI=IMAG (CKPHP ( IF , JZ) ) 

AKZPR=REAL (CKZP (IF , JZ)  ) 

AKZPI=IMAG (CKZP (IF , JZ ) ) 

AKPHNR=REAL ( CKPHN (IF, JZ) ) 

AKPHNI = IMAG ( CKPHN { I F , JZ ) ) 

AKZNR=REAL { CKZN ( IF , JZ ) ) 

AKZNI=IMAG (CKZN (IF, JZ) ) 

ALPHPR=REAL ( CLPHP ( IF , JZ ) ) 

ALPHPI = IMAG ( CLPHP ( IF , JZ ) ) 

ALZPR=REAL (CLZP (IF, JZ) ) 

ALZPI=IMAG (CLZP ( IF , JZ) ) 

ALPHNR=REAL ( CLPHN ( IF , JZ ) ) 

ALPHNI=IMAG (CLPHN (IF, JZ) ) 

ALZNR=REAL ( CLZN ( I F , JZ ) ) 

ALZNI=IMAG ( CLZN ( IF , JZ ) ) 

WRITE (31,*)  AKPHPR, ’  ‘ ,AKPHPI 

WRITE (31,*)  AKZPR, '  ' ,AKZPI 

WRITE (3 2,*)  ALPHPR, 1  ' , ALPHPI 

WRITE  ( 32  ,  *  )  ALZPR,  '  \ALZPI 
WRITE (31,*)  AKPHNR , '  ' , AKPHNI 

WRITE (31,*)  AKZNR, '  ' ,AKZNI 

WRITE { 32 , * )  ALPHNR, '  ' ,ALPHNI 

WRITE ( 32 , * )  ALZNR, 1  1 ,ALZNI 

4100  CONTINUE 

4200  CONTINUE 
CLOSE (31) 

CLOSE (32) 

END  IF 
RETURN 
END 
CC 

SUBROUTINE  INCIDNT (NA, CIWO , CIW) 

C**** ***************  ********************************* ****************** 

C  This  subroutine  sets  up  the  incident  wave  on  an  object  which  is  coated 
C  with  anisotropic  material,  or  a  perfect  conductor. 
£********************************************************************** 
INCLUDE  1  RE ALT P . INC ' 

INCLUDE  ' CMPXTP . INC ' 

INCLUDE  ' MAINDM . INC ' 

C 

DIMENSION  CIWO (KXCRT) , CIW (KCRNT) , DJNS (KNDIMl+1 ) , DJPS (KQDIMl+2 ) 
COMMON  /INPUT2 /  IE , IM, THETAI , THESINI , THECOSI , RTHEI 
COMMON  / INPUT4 /  IZ , IK, IS,NYSM 

COMMON  / INPUT3 /  RTHE , RDELT , RPHI , RDELP , NTHTAO , NPHI , THESIN , THECOS , 

+  RHPI 

COMMON  /RTHETA/  DL1COSI , DL2SINI , DLL 

COMMON  /GCONST/  DKH, DKA, HH, HA, HSQ, DASQ , HSQN, DASQN, HHSQ , HDASQ, 

+  RAH , RAHSQ , DHA , DAH 

COMMON  / CRNTDM/  NMAX , MXNG , IQMAX, IQMAX1 , IQMAX2 , IXCRNT , ICRNT , MXQEG , 


NP=NA+1 

NP1=NP+1 

NM=NA-1 

C 

IF  (IZ  .EQ.  0)  GO  TO  4000 
C 

C  This  part  computes  the  incident  field  on  an  anisotropic  object 
C 

C  Initialize  the  column  matrix  of  sum  current  on  the  anisotropic 
DO  100  IJ=1, KCRNT 
CIW(IJ) =CZERO 
100  CONTINUE 

If  the  incident  angle  is  90  degree  or  0  degree 
IF  (RTHEI  .EQ.  ZERO)  GO  TO  3000 
CALL  DBSJNS (DL2SINI ,  KNDIM1 ,  DJNS) 


DJN=DJNS  (NP) 

IF  (NA.  EQ.  0)  THEN 
DDJN=-DJNS {NP1 ) 

ELSE 

DDJN=HALF*  (DJNS  (NA)  -DJNS  (NP1 )  ) 

END  IF 

C  Use  IMSL  library  to  compute  Bessel's  function 
CALL  DBSJNS (DL1C0SI , IQMAX2 , DJPS) 

DO  600  IP=0 , KQDIM 

NIP=IP+1 

DNIP=NIP 

NIP1=NIP+1 

NIP2=NIP+2 

N11=NA+IP-1 

N12=NA+IP 

IEl=MOD (Nil , 4 ) 

IE2=MOD (N12 , 4 ) 

C 

CF1=CI2**IE1 

CF2=CI2**IE2 

IPW1=4*IP+1 

IPW2=IPW1+1 

IPW3=IPW2+1 

IPW4=IPW3+1 

C 

DJP=DJPS (NIP) 

DJP1=DJPS (NIP1 ) 

DJP2=DJPS (NIP2 ) 

IF  (IE  .EQ.  1)  THEN 
CEEPHI=CFl* (DJP+DJP2) *DDJN/DNIP 
CEHPHI=-TWO*NA*CF2*DJPl*DJN/DLL 
CEHZ=-CF2*THESINI* (DJP+DJP2) *DJN/DNIP 
ELSE 

CEEPHI=CZERO 
CEHPHI=CZERO 
CEHZ=CZERO 
END  IF 
C 

IF  (IM  .EQ.  1)  THEN 
CMEPHI=TWO*NA*CF2  *DJP1 *DJN/DLL 
CMEZ=CF2  *THESINI * (DJP+DJP2 ) *DJN/DNIP 
CMHPHI=CF1 * (DJP+DJP2) *DDJN/DNIP 
ELSE 

CHEPHI=CZERO 
CMEZ=CZERO 
CMHPHI=CZERO 
END  IF 
C 

CIW(IPWl) =TWO* (CEEPHI+*CMEPHI) 

CIW ( IPW2 ) =TWO*CMEZ 

CIW ( IPW3 ) =TWO* (CEHPHI+CMHPHI ) 

CIW ( IPW4 ) =TW0*CEHZ 
600  CONTINUE 
RETURN 
C 

3000  CALL  DBSJNS (DKH,IQMAX2, DJPS) 

DO  3200  IP=0 , KQDIM 

JP=IP+2 

DJP=DJPS ( JP) 

IE11=M0D ( IP, 4 ) 

IE12=M0D ( IP+1 , 4 ) 

CF1=CI2**IE11 

CF2=CI2**IE12 

IPW1=4*IP+1 

IPW3=IPWl+2 

C 

IF  (IE  .EQ.  1)  THEN 
CEEPHI=CF1 *DJP/DKH 


series . 
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CEHPHI=-CF2*DJP/DKH 

ELSE 

CEEPHI=CZERO 
CEHPHI=CZERO 
END  IF 
C 

IF  (IM  .EQ.  1)  THEN 
CME  PHI = CF2  *  D JP / DKH 
CMH  PH I = CF 1 *  D JP / DKH 
ELSE 

CMEPHI=CZERO 
CMHPHI=CZERO 
END  IF 
C 

CIW ( IPW1 ) =TWO* ( CEEPHI+CMEPHI ) 

CIW(IPW3)=TW0* ( CEHPHI +CMHPHI ) 

3200  CONTINUE 
RETURN 
C 

C  This  part  compute  incident  fields  on  a  perfect  conductor 
C  Initialize  the  column  matrix  of  sum  current  on  the  conductor 
4000  DO  4200  IX=1,KXCRT 
CIW0 (IX) =CZERO 
4200  CONTINUE 

If  the  incident  angle  is  90  degree  or  0  degree 
IF  (RTHEI  .EQ.  ZERO)  GO  TO  6000 
C 

CALL  DBSJNS (DL2SINI ,  MXNG+1,  DJNS) 

DJN=DJNS  (NP) 

IF  (NA  .EQ.  0)  THEN 
DDJN=-DJNS (NP1) 

ELSE 

DDJN=HALF*  (DJNS  (NA)  -DJNS  (NP1)  ) 

END  IF 
C 

CALL  DBSJNS (DL1COSI , IQMAX2 , DJPS ) 

DO  4400  IP=0 , KQDIM 

NIP=IP+1 

DNIP=NIP 

NIP1=NIP+1 

NIP2=NIP+2 

N11=NA+IP-1 

N12=NA+IP 

IEl=MOD (Nil , 4 ) 

IE2=M0D(N12, 4) 

C 

CF1=CI2**IE1 

CF2=CI2**IE2 

IPW1=2*IP+1 

IPW2=IPW1+1 

C 

DJP=DJPS (NIP) 

DJP1=DJPS (NIP1) 

DJP2=DJPS (NIP2 ) 

IF  (IE  .EQ.  1)  THEN 
CEEPHI=CF1* (DJP+DJP2) *DDJN/DNIP 
ELSE 

CEEPHI=CZERO 
END  IF 
C 

IF  (IM  .EQ.  1)  THEN 
CMEPHI=TWO*NA*CF2*DJPl*DJN/DLL 
CMEZ=CF2*THESINI* (DJP+DJP2 ) *DJN/DNIP 
ELSE 

CMEPHI=CZERO 
CMEZ=CZERO 
END  IF 
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CIWO (IPW1) =TWO* (CEEPHI+CMEPHI) 

CIWO (IPW2 ) =TWO*CMEZ 
4400  CONTINUE 
RETURN 
C 

6000  CALL  DBS JNS ( DKH , IQMAX2  ,  D JPS ) 

DO  6200  IP=0 , KQDIM 

JP=IP+2 

DJP=DJPS ( JP) 

IE11=M0D (IP, 4 ) 

IE12=MOD ( IP+1 , 4 ) 

CF1=CI2**IE11 

CF2=CI2**IE12 

IPW1=2*IP+1 

C 

IF  (IE  .EQ.  1)  THEN 
CEEPHI=CF1 *DJP/DKH 
ELSE 

CEEPHI=CZERO 
END  IF 
C 

IF  (IM  .EQ.  1)  THEN 
CMEPHI =CF2* D JP / DKH 
ELSE 

CMEPHI=CZERO 
END  IF 
C 

CIWO (IPWl ) =TWO* (CEEPHI+CMEPHI) 

6200  CONTINUE 
RETURN 
END 
C 

SUBROUTINE  XPQO 

C*  ***************************************  +  ********■*■******■************  *-*  * 
C  This  subroutine  computes  the  matrix  XN(P,Q)  for  N  =  0  following  a 
C  call  to  XPQINI.  This  matrix  is  kept  in  the  common  block: 

C  COMMON  /XPQTMP/  CXPQN 

C* ********************************************************************** 
INCLUDE  ' REALTP . INC ‘ 

INCLUDE  ' CMPXTP . INC ' 

INCLUDE  ' MAINDM . INC ' 

C 

DIMENSION  CXPQN ( KCRNT , KCRNT ) , CXRPQ ( KCRNT , KCRNT ) 

DIMENSION  CXPQN0 ( KXCRT , KXCRT ) 

DIMENSION  CGNE ( 0 : MAXPEG , 0 : MAXPEG , KNDIM1+1 ) , 

+  CGNO { 0 : MAXPOG , 0 : MAXPOG , KNDIMl + 1 ) 

DIMENSION  CDGNE ( 0 : MAXPEG , 0 : MAXPEG , KNDIMl + 1 ) , 

+  CDGNO (0 :MAXPOG, 0 :MAXPOG, KNDIMl +1) 

COMMON  /GPQTMP/  CGNE , CDGNE , CGNO , CDGNO 
COMMON  /XPQTMP/  CXPQN, CXRPQ, CXPQN0 

COMMON  / CRNTDM/  NMAX,  MXNG,  IQMAX,  IQMAXl,  IQMAX2  ,  IXCRNT,  ICRNT,  MXQEG, 

+  MXQOG 

COMMON  / GCONST /  DKH, DKA, HH, HA, HSQ, DASQ, HSQN, DASQN, HHSQ, HDASQ, 

+  RAH , RAHSQ , DHA , DAH 

COMMON  /INPUT 4/  IZ , IK, IS , NYSM 
SAVE  /XPQTMP/ , /GPQTMP/ , /GCONST/ 

C 

IF  (N  .NE.  0)  THEN 

WRITE (*,*)  ’Input  N  is  not  equal  to  0  in  XPQO.' 

WRITE { * , * )  'Execution  is  stopped.' 

STOP 
END  IF 
C 

IF  (IZ  .EQ.  0)  GO  TO  4000 
C 

C  Initialize  the  XN(P,Q)  matrix. 

C  Null  terms  will  be  skipped  later. 
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DO  200  IQ=1,KCRNT 
DO  100  IP=1 , KCRNT 
CXPQN ( IP , IQ ) =CZERO 
100  CONTINUE 

200  CONTINUE 

Form  the  XN(P,Q)  matrix  for  n  =  0  and  using  on  a  perfect  conductor 
CF3 1=DKA*CI1 

DO  1300  IQE=0 , MXQEG-1 
IQE1=IQE+1 
IQ=2*IQE 
IQX1=4*IQ+1 
IQX2=IQXl+l 
IQX3=IQX2+1 
IQX4=IQX3+1 
DQ1=IQ+1 
FQ22=DQ1/HHSQ 

DO  1100  IPE=0 , MXQEG-1 
IPE1=IPE+1 
IP=2*IPE 
IPX1=4*IP+1 
IPX2=IPX1+1 
IPX3=IPX2+1 
IPX4=IPX3+1 
DP1=IP+1 
F41=DKH/DP1 
CF11=F41*CF31 
F51=QUAR*F41*DKA 

CXPQN ( I PX1 , IQX1 ) = (CGNE ( IPE1 , IQE, 2) -CGNE ( IPE , IQE ,  2 )  ) *CF11 
CXPQN ( IPX4 , IQX1 ) = (CGNE ( IPE1 , IQE , 2 ) -CGNE ( IPE , IQE , 2 ) +HA* ( 

+  CDGNE (IPEl , IQE, 2 ) -CDGNE (IPE, IQE, 2 ) ) ) *F41 

CXPQN (IPX2, IQX2) = (TWO*FQ22*DPl*CGNO (IPE, IQE, 1) +HALF* ( 

+  CGNE (IPE, IQE1 , 1) +CGNE ( IPEl , IQE , 1 ) - 

+  CGNE (IPEl, IQE1,1) -CGNE (IPE, IQE, 1) ) ) *CF11 

CXPQN ( IPX3 , IQX2 ) = ( CDGNE ( I PE , IQE , 1 ) + CDGNE (IPEl , IQE1 , 1 ) - 
+  CDGNE ( IPEl , IQE , 1 ) -CDGNE ( IPE , IQEl , 1 ) ) *F51 

CXPQN ( IPX2 , IQX3 ) =- CXPQN ( IPX4 , IQX1 ) 

CXPQN ( IPX3 , IQX3 ) =CXPQN ( IPX1 , IQX1 ) 

CXPQN ( IPX1 , IQX4 ) = -CXPQN ( IPX3 , IQX2 ) 

CXPQN ( IPX4 , IQX4 ) =CXPQN ( IPX2 , IQX2 ) 

1100  CONTINUE 

1300  CONTINUE 

DO  2300  IQO=0 , MXQOG-1 

IQ01=IQO+l 

IQ=2*IQO+l 

IQX1=4*IQ+1 

IQX2=IQX1+1 

IQX3=IQX2+1 

IQX4=IQX3+1 

DQ1=IQ+1 

FQ22=DQ1/HHSQ 

DO  2200  IPO=0 , MXQOG-1 
IP01=IPO+l 
IP=2  *IPO+l 
IPX1=4*IP+1 
IPX2=IPX1+1 
IPX3=IPX2+1 
IPX4=IPX3+1 
DP1=IP+1 
F41=DKH/DP1 
CF11=F41*CF31 
F51=QUAR*F41 *DKA 

CXPQN ( IPX1 , IQX1 ) = ( CGNO ( IPOl , IQO , 2 ) -CGNO ( IPO , IQO , 2 ) ) *CF11 
CXPQN ( IPX4 , IQX1 ) = ( CGNO (IPOl , IQO , 2 ) -CGNO ( IPO , IQO , 2 ) +HA* ( 

+  CDGNO ( IPOl , IQO , 2 ) -CDGNO (IPO, IQO, 2 ) ) ) *F41 

CXPQN (IPX2, IQX2) = (TW0*FQ22*DP1*CGME (IPOl , IQOl , 1 ) +HALF* ( 
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+  CGNO ( IPO , IQOl , 1 ) +CGNO ( IPOl , IQO ,  1 )  - 

+  CGNO (IPOl, IQOl, 1) -CGNO (IPO, IQO, 1) ) ) *CF11 

CXPQN ( IPX3 , IQX2 )  =  ( CDGNO ( IPO , IQO , 1 ) + CDGNO ( IPOl , IQOl ,  1 )  - 
+  CDGNO (IPOl, IQO, 1) -CDGNO (IPO, IQOl, 1) ) *F51 

CXPQN ( IPX2 , IQX3 ) = -CXPQN ( IPX4 , IQXl ) 

CXPQN ( IPX3 , IQX3 ) =CXPQN ( IPX1 , IQXl ) 

CXPQN ( IPX1 , IQX4 ) = -CXPQN ( IPX3 , IQX2 ) 

CXPQN (IPX4 , IQX4 ) =CXPQN (IPX2 , IQX2 ) 

2200  CONTINUE 

2300  CONTINUE 
RETURN 

C  Initialize  the  XN(P,Q)  matrix. 

C  Null  terms  will  be  skipped  later. 

C 

4000  DO  4200  IQ=1,KXCRT 

DO  4100  IP=1 , KXCRT 
CXPQN 0 ( IP , IQ) =CZERO 
4100  CONTINUE 

4200  CONTINUE 

Form  the  XN(P,Q)  matrix  for  n  =  0  and  using  on  an  anisotropic  coat. 
CF3 1=DKA*CI1 
C 

DO  4500  IQE=0 , MXQEG-1 

IQE1=IQE+1 

IQ=2*IQE 

IQX1=2*IQ+1 

IQX2=IQX1+1 

DQ1=IQ+1 

FQ22=DQ1/HHSQ 

DO  4400  IPE=0, MXQEG-1 
IPE1=IPE+1 
IP=2*IPE 
IPX1=2*IP+1 
IPX2=IPX1+1 
DP1=IP+1 
F41=DKH/DP1 
CF11=F41*CF31 
F51=QUAR*F41*DKA 

CXPQN0  ( IPX1 ,  IQXl )  =  ( CGNE  ( IPEl ,  IQE ,  2  )  -CGNE  (IPE,  IQE,  2  )  )  *CF11 
CXPQN0 ( IPX2 , IQX2 ) = (TWO*FQ22*DPl*CGNO (IPE, IQE, 1) +HALF* ( 

+  CGNE ( IPE , IQE1 , 1 ) +CGNE ( IPEl , IQE , 1 ) - 

+  CGNE ( IPEl, IQE1,1) -CGNE (IPE, IQE, 1) ) ) *CF11 

4400  CONTINUE 

4500  CONTINUE 

DO  5300  IQO=0 , MXQOG-1 

IQ01=IQO+l 

IQ=2*IQO+l 

IQX1=2*IQ+1 

IQX2=IQX1+1 

DQ1=IQ+1 

FQ22=DQ1/HHSQ 

DO  5200  IPO=0, MXQOG-1 
IP01=IPO+l 
IP=2*IPO+l 
IPX1=2*IP+1 
IPX2=IPX1+1 
DP1=IP+1 
F41=DKH/DP1 
CF11=F41*CF3 1 
F51=QUAR*F41*DKA 

CXPQN 0 ( I PX1 , IQXl )  =  ( CGNO ( IPOl , IQO , 2 ) -CGNO { IPO ,  IQO ,  2  )  )  *CF11 
CXPQNO ( IPX2 , IQX2 ) = (TW0*FQ22*DP1*CGNE (IPOl , IQOl , 1 ) +HALF* ( 

+  CGNO ( IPO , IQOl , 1 ) +CGNO (IPOl , IQO , 1) - 

+  CGNO ( I POl , IQOl , 1 ) -CGNO ( IPO , IQO , 1 ) ) ) *CF11 

5200  CONTINUE 

5300  CONTINUE 
RETURN 
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END 

C 

SUBROUTINE  XPQN(NIN) 

Q* ********************************************************************** 

C  This  subroutine  calls  the  subroutine  GDN  to  update  G(P,Q,N)  for 
C  N  =  NIN+1,  then  forms  the  matrix  XN(P,Q)  for  N  =  NIN  >  0.  This  matrix 
C  is  kept  in  the  common  block: 

C  COMMON  /XPQTMP/  CXPQN 

C* ********************************************************************** 

C  WARNING: 

C  IT  IS  ASSUMED  THAT  XPQINI  AND  XPQN  FOR  N  FROM  1  TO  NIN-1  HAVE  BEEN 
C  CALLED  SO  THAT  G(P,Q,N)  AND  ITS  DERIVATIVE  FOR  N=NIN-1  AND  N=NIN 
C  HAVE  BEEN  STORED  IN  THE  COMMON  BLOCK  /GPQTMP/  WITH  PROPER  N-INDICES. 

INCLUDE  1  RE ALT P . INC  1 
INCLUDE  1 CMPXTP . INC  1 
INCLUDE  ' MAINDM . INC  1 
C 

DIMENSION  CXPQN ( KCRNT ,  KCRNT ) , CXRPQ (KCRNT , KCRNT ) 

DIMENSION  CXPQNO (KXCRT, KXCRT) 

DIMENSION  CGNE ( 0 : MAXPEG , 0 : MAXPEG , KNDIM1 + 1 ) , 

+  CGNO ( 0 : MAXPOG , 0 : MAXPOG , KNDIM1+1 ) 

DIMENS I ON  CDGNE ( 0 : MAXPEG , 0 : MAXPEG , KNDIM1 +1 ) , 

+  CDGNO ( 0 : MAXPOG , 0 : MAXPOG / KNDIMl+1 ) 

COMMON  /GPQTMP/  CGNE , CDGNE , CGNO , CDGNO 
COMMON  /XPQTMP/  CXPQN, CXRPQ, CXPQNO 

COMMON  / CRNTDM/  NMAX, MXNG, IQMAX, IQMAXl , IQMAX2 , IXCRNT, ICRNT, MXQEG, 

+  MXQOG 

COMMON  / GCONST/  DKH , DKA , HH , HA , HSQ , DASQ , HSQN , DASQN , HHSQ , HDASQ , 

+  RAH , RAHSQ , DHA , DAH 

COMMON  /NCONST/  DN , DNH , N 
COMMON  /INPUT 4 /  IZ , IK, IS , NYSM 
SAVE  /XPQTMP/ , /GPQTMP/ , /GCONST/ 

C 

N=NIN 
NA=ABS (N) 

IF  {NA  . LT .  1)  THEN 

WRITE ( * , * )  'Input  ABS (N)  is  less  than  1  in  XPQN.1 
WRITE ( * , * )  'Execution  is  stopped.' 

STOP 
END  IF 
C 

N0I=WA+1 

NPI=N0I+1 

NMI=NA 

C 

DN0=NA 

DNP=N0I 

DNM=NA-1 

C 

IF  (IZ  .EQ.  0)  GO  TO  4000 
C 

C  Initialize  the  XN(P,Q)  matrix. 

C  Null  terms  will  be  skipped  later. 

C . 

DO  800  IQ=1, KCRNT 
DO  700  IP=1, KCRNT 
CXPQN (IP, IQ) =CZERO 
700  CONTINUE 

800  CONTINUE 

Form  the  XN(P,Q)  matrix. 

CF3 1=DKA*CI1 
F21=TWO*DNO 
FN1 1 =F2 1 * DN 0 / DASQ 
C 

DO  1300  IQE=0 , MXQEG-1 
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IQE1=IQE+1 

IQ=2*IQE 

IQXl=4*IQ+l 

IQX2=IQX1+1 

IQX3=IQX2+1 

IQX4=IQX3+1 

DQ1=IQ+1 

FQ12=DN0*DQ1 

FQ22=DQ1/HHSQ 

DO  1100  IPE=0 , MXQEG-1 
IPE1=IPE+1 
IP=2  *IPE 
IPX1=4*IP+1 
IPX2=IPXl+l 
IPX3=IPX2+1 
IPX4=IPX3+1 
DP1=IP+1 
F41=HH/DP1 
CF11=F41*CF31 
F51=HALF*DKA*F41 

CXPQN ( IPX1 , IQX1 )  =  ( (CGNE ( IPE , IQE ,  NOI ) -CGNE ( IPEl , IQE f  NOI ) ) *FN11+ 
+  CGNE { IPEl , IQE , NMI ) -CGNE (IPE , IQE, NMI ) + 

+  CGNE (IPEl, IQE, NPI) -CGNE (IPE, IQE,NPI) ) *CF11 

CXPQN ( IPX4 , IQXl ) = ( ( CGNE ( I PE , IQE , NMI ) -CGNE ( IPEl , IQE , NMI ) ) *DNM+ 

+  { CGNE ( IPEl , IQE , NPI ) -CGNE ( IPE , IQE , NPI ) ) *DNP+ 

+  HA* ( CDGNE ( IPEl , IQE , NMI ) -CDGNE ( IPE , IQE , NMI ) + 

+  CDGNE (IPEl, IQE, NPI) -CDGNE (IPE, IQE, NPI) ) ) *F41 

CXPQN ( IPX2 , IQX2 )  =  (FQ22  *DP1 *CGNO ( IPE , IQE , NOI )  + 

+  CGNE ( IPE , IQE1 , NO I ) +CGNE ( IPEl , IQE , NO I ) - 

+  CGNE (IPEl , IQE1 , NOI ) -CGNE (IPE,IQE,N0I))*CF11 

CXPQN ( IPX3 , IQX2 ) = ( CDGNE ( IPE , IQE , NO I ) + CDGNE ( IPEl , IQE1 , NOI ) - 
+  CDGNE (IPEl, IQE, NO I) -CDGNE ( IPE , IQE1 , NOI ) ) *F51 

CXPQN ( IPX2 , IQX3 ) =-CXPQN ( IPX4 , IQXl ) 

CXPQN ( IPX3 , IQX3 ) =CXPQN ( IPX1 , IQXl ) 

CXPQN ( IPX1 , IQX4 ) = -CXPQN ( IPX3 , IQX2 ) 

CXPQN ( IPX4 , IQX4 ) =CXPQN ( IPX2 , IQX2 ) 

1100  CONTINUE 

DO  1200  IPO=0 , MXQOG-1 
IP01=IP0+1 
IP=2*IP0+1 
IPX1=4*IP+1 
IPX2=IPXl+l 
IPX3=IPX2+1 
IPX4=IPX3+1 
DP1=IP+1 
F12=FQ12/DP1 

CXPQN  ( IPX2  ,  IQXl )  =-F21*CGNE  { IP01 ,  IQE ,  NOI ) 

CXPQN ( IPX3 , IQXl ) = (CGNE ( IPOl , IQE , NMI ) -CGNE ( IPOl , IQE , NPI ) ) *CF3 1 
CXPQN ( IPX1 , IQX2 ) = (CGNO ( IPOl , IQE , NOI ) -CGNO ( IPO , IQE , NOI ) ) *F12 
CXPQN { IPX1 , IQX3 ) = -CXPQN ( IPX3 , IQXl ) 

CXPQN { IPX4, IQX3)=CXPQN(IPX2, IQXl) 

CXPQN ( IPX3 , IQX4 ) =CXPQN ( IPX1 , IQX2 ) 

1200  CONTINUE 

1300  CONTINUE 

DO  2300  IQO=0, MXQOG-1 

IQ01=IQ0+1 

IQ=2*IQO+l 

IQX1=4*IQ+1 

IQX2=IQX1+1 

IQX3=IQX2+1 

IQX4=IQX3+1 

DQ1=IQ+1 

FQ12=DN0*DQ1 

FQ22=DQ1/HHSQ 

DO  2100  IPE=0 , MXQEG-1 
IPE1=IPE+1 
IP=2*IPE 
IPX1=4*IP+1 
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IPX2=IPXl+l 
IPX3=IPX2+1 
IPX4=IPX3+1 
DP1=IP+1 
F12=FQ12 /DPI 

CXPQN { IPX2 , IQX1 ) =F21* CGNO { I PE , IQO ,  NO I ) 

CXPQN ( IPX3 , IQX1)= {CGNO (IPE, IQO, NMI) -CGNO (I PE, IQO, NPI) ) *CF31 
CXPQN (IPX1 , IQX2 ) = (CGNE (IPE1 , IQOl , NOI) -CGNE (IPE, IQOl,N0l) ) *F12 
CXPQN { IPX1 , IQX3 ) = -CXPQN ( IPX3 , IQX1 ) 

CXPQN ( IPX4 , IQX3 ) =CXPQN ( IPX2 , IQX1 ) 

CXPQN { IPX3 , IQX4 ) =CXPQN ( IPX1 , IQX2 ) 

2100  CONTINUE 

DO  2200  IPO=0 , MXQOG-1 
IP01=IP0+1 
IP=2*IPO+l 
IPX1=4*IP+1 
IPX2=IPX1+1 
IPX3=IPX2+1 
IPX4=IPX3+1 
DP1=IP+1 
F41=HH/DP1 
CF11=F41*CF31 
F51=HALF*DKA*F41 

CXPQN ( IPX1 , IQX1 ) = ( ( CGNO ( IPO, IQO, NO I ) -CGNO (IPOl, IQO, NOI) )*FN11+ 
+  CGNO { IPOl , IQO , NMI ) -CGNO ( IPO , IQO , NMI ) + 

+  CGNO {IPOl, IQO, NPI) -CGNO (IPO, IQO,NPI) ) *CF11 

CXPQN { IPX4 , IQX1 ) = ( (CGNO { IPO , IQO , NMI ) -CGNO ( IPOl , IQO , NMI ) ) *DNM+ 

+  {CGNO {IPOl, IQO, NPI) -CGNO (IPO, IQO, NPI) ) *DNP+ 

+  HA* ( CDGNO ( IPOl , IQO , NMI ) -CDGNO ( IPO , IQO , NMI ) + 

+  CDGNO (IPOl, IQO, NPI) -CDGNO (IPO, IQO, NPI) ) ) *F41 

CXPQN { IPX2 , IQX2 ) = (FQ22*DP1*CGNE { IPOl , IQOl , NOI ) + 

+  CGNO ( IPO , IQOl , NOI ) +CGNO ( IPOl , IQO , NOI ) - 

+  CGNO (IPOl, IQOl, NOI) -CGNO (IPO, IQO, NOI) ) *CF11 

CXPQN { IPX3 , IQX2 ) = ( CDGNO ( IPO , IQO , NO I ) +CDGNO ( IPOl , IQOl , NO I ) - 
+  CDGNO (IPOl , IQO , NOI ) -CDGNO ( IPO , IQOl , NO I ) ) *F51 

CXPQN (IPX2 , IQX3 ) = -CXPQN ( IPX4 , IQXl ) 

CXPQN ( IPX3 , IQX3 ) =CXPQM ( IPX1 , IQXl ) 

CXPQN ( IPX1 , IQX4 ) = -CXPQN ( IPX3 , IQX2 ) 

CXPQN { IPX4 , IQX4 ) =CXPQN ( IPX2 , IQX2 ) 

2200  CONTINUE 

2300  CONTINUE 
RETURN 

C  Initialize  the  XN(P,Q)  matrix. 

C  Null  terms  will  be  skipped  later. 

C 

4000  DO  4800  IQ=1,KXCRT 

DO  4700  IP=1 , KXCRT 
CXPQN0 (IP, IQ) =CZERO 
4700  CONTINUE 

4800  CONTINUE 
C 

C  Form  the  XN(P,Q)  matrix. 

CF3 l=DKA*CIl 
F21=TWO*DNO 
FN11=F21*DN0/DASQ 
C 

DO  5300  IQE=0 , MXQEG-1 

IQE1=IQE+1 

IQ=2*IQE 

IQX1=2*IQ+1 

IQX2=IQX1+1 

DQ1=IQ+1 

FQ12=DN0*DQ1 

FQ22=DQ1/HHSQ 

DO  5100  IPE=0, MXQEG-1 
IPE1=IPE+1 
IP=2*IPE 
IPX1=2*IP+1 
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IPX2=IPX1+1 

DP1=IP+1 

F41=HH/DP1 

CFll=F41*CF31 

F51=HALF*DKA*F41 

CXPQNO  (IPX1,  IQX1)  =  (  ( CGNE  ( IPE ,  IQE ,  NO  I )  -CGNE  ( IPEl ,  IQE ,  NOI )  )  *FN11  + 

+  CGNE ( I PEI , IQE , NMI ) -CGNE (IPE, IQE, NMI ) + 

+  CGNE ( IPEl , IQE , NPI ) -CGNE ( IPE , IQE , NPI ) ) *CFll 

CXPQNO ( IPX2 , IQX2 ) = (FQ22*DP1*CGN0 (IPE, IQE, NOI ) + 

+  CGNE ( IPE , IQE1 , NO I ) +CGNE (IPEl, IQE , NO I ) - 

+  CGNE (IPEl, IQE1, NOI) -CGNE (IPE, IQE, NOI) ) *CF11 

5100  CONTINUE 

DO  5200  IPO=0 , MXQOG-1 
IP01=IP0+1 
IP=2*IP0+1 
IPX1=2*IP+1 
IPX2=IPX1+1 
DP1=IP+1 
F12=FQ12/DP1 

CXPQNO ( IPX2 , IQX1 ) =F21*CGNE ( IP01 , IQE, NOI ) 

CXPQNO ( IPX1 , IQX2 ) = ( CGNO ( IP01 , IQE , NOI ) -CGNO ( IPO , IQE , NOI ) ) *F12 
5200  CONTINUE 

5300  CONTINUE 

DO  6300  IQO=0, MXQOG-1 

IQ01=IQ0+1 

IQ=2*IQO+l 

IQX1=2*IQ+1 

IQX2=IQX1+1 

DQ1=IQ+1 

FQ12=DN0*DQl 

FQ22=DQ1/HHSQ 

DO  6100  IPE=0 , MXQEG-1 
IPE1=IPE+1 
IP=2*IPE 
IPX1=2*IP+1 
IPX2=IPXl+l 
DP1=IP+1 
F12=FQ12 /DPI 

CXPQNO (IPX2, IQX1)=F21*CGN0 (IPE, IQO,NOI) 

CXPQNO ( IPX1 , IQX2 ) = ( CGNE ( IPEl , IQOl , NOI ) -CGNE ( IPE , IQOl , NO I ) ) *F12 
6100  CONTINUE 

DO  6200  IPO=0, MXQOG-1 
IP01=IP0+1 
IP=2*IP0+1 
IPX1=2*IP+1 
IPX2=IPX1+1 
DP1=IP+1 
F41=HH/DP1 
CF11=F41*CF3 1 
F51=HALF^DKA*F41 

CXPQNO ( I PX1 , I QXl ) = ( ( CGNO ( I PO , I QO , NO I ) -CGNO ( I POl , IQO , NO I ) ) *FN1 1 + 

+  CGNO ( IPOl , IQO , NMI ) -CGNO ( IPO , IQO , NMI ) + 

+  CGNO (IPOl, IQO, NPI) -CGNO (IPO, IQO, NPI) ) *CF11 

CXPQNO (IPX2 , IQX2 ) = (FQ22*DP1*CGNE (IPOl , IQOl, NOI ) + 

+  CGNO (IPO, IQOl, NOI) +CGN0(IP01, IQO, NOI) - 

+  CGNO (IPOl, IQOl, NOI) -CGNO (IPO, IQO, NOI) ) *CF11 

6200  CONTINUE 

6300  CONTINUE 
RETURN 
END 
C 

SUBROUTINE  ESCFAR (NIN, CKLNO , CKLN, CETHN, CEPHN) 

Q***********************************************************************^******* 

C  This  subroutine  computes  the  scattered  fields  in  far  zone 
C 

INCLUDE  ' REALTP . INC ' 

INCLUDE  ' CMPXTP . INC ’ 

INCLUDE  ' MAINDM . INC ' 
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c 

DIMENSION  CKLNO (KXCRT) 

DIMENSION  CKLN(KCRNT) , DJNS (KNDIM1 ) , DJPS (KQDIM1+2 ) 

COMMON  / CRNTDM/  NMAX , MXNG , IQMAX , IQMAX1 , IQMAX2 , IXCRNT , ICRNT , MXQEG , 
+  MXQOG 

COMMON  / GCONST/  DKH , DKA, HH, HA, HSQ, DASQ , HSQN, DASQN, HHSQ, HDASQ , 

+  RAH , RAH  SQ , DHA , DAH 

COMMON  / INPUT 3 /  RTHE , RDELT , RPHI , RDELP , NTHTAO , NPHI , THESIN , THECOS , 

+  RHPI 

COMMON  / INPUT4 /  IZ, IK, IS,NYSM 
C 

N=NIN 

DN1=N 

NA=ABS(NIN) 

NP=NA+1 
NPl=NP+l 
PHI1=RPHI*N 
CI2N=CI2*  *N 

CENP=CONE*COS < PHI1 ) +CIl*SIN ( PHI1 ) 

CETHN=CZERO 

CEPHN=CZERO 

C 

IF  ((RTHE  .EQ.  0.)  .OR.  (RTHE  .E Q.  PI))  GO  TO  2000 
C 

TCOT=THECOS / THES IN 
DL2SIN=DKA*THESIN 
DLlCOS=DKH*THECOS 
DNL2 =DN1 *  TCOT / DKA 

CALL  DBSJNS (DL2SIN,  KNDIMl ,  DJNS) 

C 

DJN=DJNS (NP) 

IF  (N  .EQ.  0)  THEN 
DDJN=-DJNS (NPl ) 

ELSE 

DDJN=HALF*  (DJNS  (NA)  -DJNS  (NPl)  ) 

ENDIF 

C 

KN=NA/2 

KNP=NP/2 

IF  { (N  .LT.  0)  .AND.  (KNP  .GT.  KN)  )  THEN 
DJN=-DJN 
DDJN=-DDJN 
ENDIF 
C 

CALL  DBSJNS (DL1COS,  KQDIM1+2,  DJPS) 

C 

IF  (IZ  .EQ.  0)  THEN 
DO  100  IP=0, IQMAX 
INP=IP+1 
NP2=INP+2 
N11=IP 
Nl2=IP+2 
IEl=MOD (Nil , 4 ) 

CF1=CI2**IE1 
IPW1=2*IP+1 
IPW2=IPW1+1 
DJP=DJPS ( INP ) 

DJP2=DJPS(NP2) 

DJPH=HALF* (DJP+DJP2 ) 

CETHN=CETHN+CI2N*DHA*CENP*CI1*CF1*DJN* (DNL2* CKLNO (IPW1) *DJP 
+  -THESIN* CKLNO (IPW2) *DJPH) 

CEPHN=CEPHN-CI2N*DHA*CENP*DDJN*CKLN0 (IPW1) *CF1*DJP 
100  CONTINUE 
ELSE 

DO  200  IP=0, IQMAX 
INP=IP+1 
NP2=INP+2 
N11=IP 
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Nl2=IP+2 
IE1=M0D (Nil , 4 ) 

CFl=CI2  **IE1 

IPW1=4*IP+1 

IPW2=IPW1+1 

IPW3=IPW2+1 

IPW4=IPW3+1 

DJP=DJPS (INP) 

DJP2=DJPS (NP2 ) 

DJPH=HALF* (DJP+DJP2 ) 

CETHN=CETHN+CI2N*DHA*CENP* ( (CI1*DJN*DNL2 *CKLN (IPW1) -DDJN 
+  *  CKLN ( I PW3 ) ) *DJP-CI1*DJN*THESIN*CKLN ( IPW2 ) *DJPH) *CF1 

CEPHN=CEPHN-CI2N*DHA*CENP* ( (DDJN* CKLN ( IPWl ) +CI1*DJN*DNL2 
+  *CKLN(IPW3) ) *DJP-CI1*DJN*THESIN*CKLN ( IPW4 ) *DJPH) *CF1 

200  CONTINUE 
END  IF 
RETURN 
C 

2000  IF  (NA  .NE.  1)  THEN 

GO  TO  3000 
END  IF 

CALL  DBS JNS ( DKH , KQDIM1 , D JPS ) 

C 

IF  (IZ  .EQ.  0)  THEN 
DO  2100  IP=0 / IQMAX 
INP=IP+1 
IPW1=2*IP+1 
JP=MOD ( IP , 4 ) 

CF1=CI1* * JP 
CF2=CI2  ** JP 
C 

DJP=DJPS (INP) 

IF  (RTHE  .EQ.  0.)  THEN 
IF  (N  .LT.  0)  THEN 
CETHN=CETHN-DAH*DJP*CF2  *CKLN0 ( IPWl ) 
CEPHN=CEPHN+DAH*DJP*CF2*CI1*CKLN0 ( IPWl ) 

ELSE 

CETHN=CETHN+DAH*DJP*CF2  *CKLN0 (IPWl ) 

CEPHN=CEPHN+DAH*DJP*CF2  *CI1*CKLN0 ( IPWl ) 

END  IF 
END  IF 

IF  (RTHE  .EQ.  PI)  THEN 
IF  (N  .LT.  0)  THEN 
CETHN=CETHN-DAH  *  D JP  *  CF 1 *  CKLNO ( IPWl ) 

CEPHN= CE  PHN+DAH  *  D JP  *CF1*CI1*  CKLNO ( I PW1 ) 

ELSE 

CETHN=CETHN+DAH*DJP*CF1*CKLN0 ( IPWl ) 
CEPHN=CEPHN+DAH*DJP*CF1*CI1* CKLNO ( IPWl ) 

END  IF 
END  IF 

2100  CONTINUE 
ELSE 

DO  2200  IP=0, IQMAX 
INP=IP+1 
IPW1=4*IP+1 
IPW3=IPWl+2 
JP=MOD ( IP , 4 ) 

CF1=CI1**JP 

CF2=CI2**JP 

C 

DJP=DJPS ( INP) 

IF  (RTHE  .EQ.  0.)  THEN 
IF  (N  .LT.  0)  THEN 

CETHN=CETHN-DAH*DJP*CF2  * (CKLN ( IPWl ) -CI1*CKLN (IPW3 ) ) 

CEPHN= CEPHN+DAH  *DJP*CF2*(CIl*  CKLN ( I PW1 ) +CKLN ( I PW3 ) ) 

ELSE 

CETHN=CETHN+DAH*DJP*CF2  * (CKLN ( IPWl ) +CI1* CKLN ( I PW3 ) ) 
CEPHN=CEPHN+DAH*DJP*CF2* (CI1 * CKLN ( IPWl ) -CKLN (I PW3 ) ) 
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END  IF 
END  IF 

IF  (RTHE  .EQ.  PI)  THEN 
IF  (N  .LT.  0)  THEN 

CETHN=CETHN-DAH*DJP*CF1* (CKLN (IPWl ) +CI1*CKLN (IPW3 ) ) 
CEPHN=CEPHN+DAH*DJP*CF1* (CI1*CKLN ( IPWl ) -CKLN (IPW3 ) ) 

ELSE 

CETHN=CETHN+DAH*DJP*CF1* (CKLN ( IPWl ) -CI1*CKLN ( IPW3 ) ) 
CEPHN=CEPHN+DAH*DJP*CF1* (CI1*CKLN ( IPWl ) +CKLN ( IPW3 ) ) 

END  IF 
END  IF 

2200  CONTINUE 
END  IF 
3000  RETURN 
END 
C 

SUBROUTINE  KLCRNT (DN, CKLN, CIW) 

C*  ************************************  *  ***********************************  * 
C  This  subroutine  computes  equivalent  currents  K  and  L  on  inner  and  outer 
C  surfaces  respectively. 

INCLUDE  ' REALTP . INC ' 

INCLUDE  ' CMPXTP . INC  * 

INCLUDE  ' MAINDM . INC  1 
C 

DIMENSION  CKLN ( KCRNT ), CIW (KCRNT) 

DIMENSION  CKLNP (KCRNT) , CKLNN ( KCRNT) 

DIMENS I ON  CRPQ3 1 ( KCRNT , KCRNT ) , CRPQ3 2 ( KCRNT , KCRNT ) 

DIMENSION  CKPHP ( 61 , 61) ,CKZP(61, 61) , CLPHP ( 61 , 61 ) , 

+  CLZP(61, 61) ,  CKPHN (61,61) ,CKZN(61, 61) , 

+  CLPHN ( 61 , 61) ,CLZN(61, 61) 

COMMON  /CKLMTX/  CKPHP , CKZP/ CLPHP, CLZP, CKPHN, CKZN, CLPHN, CLZN 
COMMON  /XPQTMP2/  CRPQ3 1 , CRPQ32 

COMMON  / CRNTDM/  NMAX , MXNG , IQMAX , IQMAX1 , IQMAX2 , IXCRNT , ICRNT , MXQEG , 

+  MXQOG 

COMMON  / INPUT3 /  RTHE , RDELT , RPHI , RDELP , NTHTAO , NPHI , THESIN , THECOS , 

+  RHP  I 

COMMON  /INPUT4 /  IZ , IK , IS , NYSM 
C 

DO  400  IX=1 , ICRNT 
CKLNP ( IX) =CZERO 
CKLNN (IX) =CZERO 
400  CONTINUE 

PHI 1=DN* RPHI 

CPHI=CONE*COS ( PHI1 ) +CI1*SIN (PHI1 ) 

C 

DO  600  IX=1, ICRNT 
DO  500  IY=1, ICRNT 

CKLNP (IX) =CKLNP (IX) +CRPQ31 (IX, IY) *CKLN(IY) 

500  CONTINUE 

600  CONTINUE 
C 

DO  800  IX=1, ICRNT 
DO  700  IY=1, ICRNT 

CKLNN (IX) = CKLNN ( IX) +CRPQ32 (IX, IY) *CKLN (IY) 

700  CONTINUE 

800  CONTINUE 
C 

DO  1600  IF=-3 0 , 3  0 
DF=IF 

DPHI=DF*PI/32 . DO 
DNPHI=DN*DPHI 

CPHI=CONE‘*'DCOS  (DNPHI)  +CI1*DSIN  (DNPHI ) 

DO  1500  JZ=-30 , 30 
DJZ=JZ 

DZ-DJZ/32 . DO 
DV=DACOS (DZ) 

SINV=DSIN(DV) 
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c 

CKPHP ( IF , JZ ) =CZERO 
CKZP ( IF , JZ ) =CZERO 
CLPHP(IF, JZ)=CZERO 
CLZP ( IF / JZ ) =CZERO 
CKPHN(IF, JZ) =CZERO 
CKZN ( IF , JZ ) =CZERO 
CLPHN ( IF , JZ ) =CZERO 
CLZN ( IF / JZ ) =CZERO 
C 

DO  1400  IP=0 , IQMAX 

P=IP 

P1=IP+1 

DPV=P*DV 

DP1V=P1*DV 

DCSPV=DCOS (DPV) 

SINP1V=DSIN (DP1V) 

IPX1=IP*4+1 

IPX2=IPX1+1 

IPX3=IPX2+1 

IPX4=IPX3+1 

CKPHP ( IF , JZ ) = CKPHP ( IF , JZ ) +CPHI * ( CKLNP { I PX1 ) +CIW { I PX4 ) ) 

+  *DCSPV/PI/ SINV 

CKZP (IF, JZ) =CKZP (IF, JZ) +CPHI* { CKLNP { IPX2 ) -CIW { IPX3 ) ) *SINP1V/PI 
CLPHP { IF , JZ ) =CLPHP { IF , JZ ) +CPHI * { CKLNP { IPX3 ) -CIW ( IPX2 ) ) 

+  *DCSPV/PI/ SINV 

CLZP (IF, JZ) =CLZP ( IF , JZ ) +CPHI * { CKLNP { I PX4 ) +CIW ( IPX1 ) ) *  SINP1V/ PI 
CKPHN { IF , JZ ) =CKPHN ( IF , JZ ) +CPHI * (CKLNN(IPXl) -CIW(IPX4) ) 

+  *DCSPV/PI/SINV 

CKZN ( IF , JZ ) =CKZN ( IF , JZ ) +CPHI * ( CKLNN ( IPX2 ) +CIW ( IPX3 ) ) *  SINP1V/ PI 
CLPHN ( IF , JZ ) =CLPHN ( IF , JZ ) +CPHI * ( CKLNN ( IPX3 ) +CIW ( IPX2 ) ) 

+  *DCSPV/PI/SINV 

CLZN (IF, JZ) =CLZN (IF, JZ) +CPHI* (CKLNN (IPX4 ) -CIW { IPX1 ) ) *SINP1V/PI 
1400  CONTINUE 

1500  CONTINUE 

1600  CONTINUE 
RETURN 
END 


J.  SUBROUTINE  RCSPAREA 


SUBROUTINE  RCSPAREA (CESTH , CESPH ) 

c**************************************************************************** 
C  This  subroutine  computes  cross  section  per  projected  area  in  all  direction 
C 

INCLUDE  ' REALTP . INC ’ 

INCLUDE  ' CMPXTP . INC ' 

C 

REAL*  4  ARCSPPA1 , ARCSPPA2 , APHASE1 , APHASE2 

COMMON  / GCONST/  DKH , DKA, HH , HA, HSQ, DASQ, HSQN, DASQN, HHSQ, HDASQ, 

+  RAH , RAHSQ , DHA , DAH 

COMMON  / INPUT3 /  RTHE, RDELT , RPHI , RDELP , NTHTAO, NPHI , THESIN, THECOS , 

+  RHP  I 

C 

ESTH=ABS (CESTH) 

ESPH=ABS (CESPH) 

ESTHSQ=ESTH*ESTH 

ESPHSQ=ESPH*ESPH 

IF  (RTHE  .EQ.  RHPI )  THEN 
RCSPPAl=PI*ESTHSQ/DKH/DKA 
RCSPPA2=PI*ESPHSQ/DKH/DKA 

ELSE  IF  ((RTHE  .EQ.  ZERO)  .OR.  (RTHE  . EQ .  PI))  THEN 
RCSPPA1  =  4 . DO  *ESTHSQ/DASQ 
RCSPPA2=4 . D0*ESPHSQ/DASQ 
ELSE 


ill 


AP=ABS (QUAR*DASQ*THECOS) +ABS <DKH*DKA*THESIN/PI ) 
RCSPPA1=ESTHSQ/AP 
RCSPPA2=ESPHSQ/AP 
END  IF 

ARC  S  P  PAl = RC S  P PA1 
ARCSPPA2=RCSPPA2 
C 

ER1=REAL { CESTH ) 

EI1=IMAG (CESTH) 

ER2  =REAL ( CESPH ) 

EI2=IMAG (CESPH) 

PHASEl =DATAN2 (Ell , ERl ) 

PHASE2  =DATAN2 (EI2 , ER2 ) 

APHASE1= PHASEl 
APHASE2=PHASE2 

WRITE  (21,*)  ARCSPPAl , APHASE1 , ARCSPPA2 , APHASE2 
C 

RETURN 

END 
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